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1.  Introduction 

Ray  theory,  modified  ray  theory,  and  normal  mode  theory 
are  three  different  approaches  for  calculating  acoustic  intensity 
at  a  given  point  in  the  ocean  once  the  source  location  and 
prevailing  sound  velocity  profile  have  been  specified.  Each 
has  its  own  limitations  and  regions  of  intended  use,  Ray  theory 
(1-4)  has  traditionally  been  the  ."ost  widely  used  approach 
for  propagation  loss  calculations  for  both  cw  (5-?)  and  explosive 
sources  (8-11),  Its  popularity  rests  on  both  its  relative  ease 
of  calculation  and  its  intuitively  appealing  description  of  the 
channeling  of  energy  along  ray  paths  and  within  ray  tubes  of 
varying  cross-section.  More  recently,  modified  ray  theory  (12-17) 
has  seen  increased  use  as  a  means  of  enhancing  ray  theory  results 
by  yielding  predictions  on  caustics  and  in  the  adjacent  caustic- 
related  shadow  zones,  regions  in  which  ray  theory  is  invalid. 
Normal  mode  theory  (18-23),  actually  the  mos  general  of  the 
three,  has  suffered  until  recently  from  a  lack  of  the  large, 
high  speed  computers  needed  for  the  solution  of  many  realistic 
propagation  problems.  It  has  also  .suffered  from,  and  continues 
to  suffer  from,  the  absence  of  a  clear,  intuitive  description 
of  the  normal  modes  that  could  compete  with  ray  diagrams 
(Figure  1.1)  in  explaining  how  energy  gets  from  source  to 
receiver.  With  the  advent  of  larger  computers,  however,  normal 
mode  theory  has  seen  development  along  several  lines  of  approach, 
and  it  is  fast  becoming  a  widely  used  technique  with  as  much 
of  a  following  as  ray  theory. 


1 


NOLTR  7^-95 


As  this  continuing  development  has  progressed,  more  effort 
has  gone  into  comparisons  of  normal  mode  theory  with  ray  theory, 
modified  ray  theory,  and  completely  different  approaches  (24). 

The  object  has  been  to  determine  the  conditions  under  which 
each  theory  is  the  appropriate  choice  for  intensity  calculations. 
From  relatively  small  scale  efforts  at  comparing  two  theories 
with  each  other  (25-2?),  efforts  have  progressed  to  complete 
detailed  comparisons  of  mode,  ray,  and  modified  ray  theories 
for  an  analytic  sound  velocity  profile  (16).  In  this  recent 
article,  Pedersen  and  Gordon  demonstrated  how  normal  mode  theory 
could  be  used  to  determine  the  validity  of  ray  theory  and 
modified  ray  theory  in  the  region  where  all  three  overlap,  near 
a  caustic.  They  did  this  comparison  for  an  analytic,  monotonically 
decreasing  sound  velocity  profile  yielding  a  close-in  caustic, 
and  they  suggested  that  it  would  be  useful  to  do  the  same  comparison 
for  a  deep  ocean  profile  yielding  a  convergence  zone  caustic. 

For  some  time, this  author  has  been  involved  in  studies  of 
convei-gence  zone  caustics  and  their  effect  on  propagating 
shock  waves  from  underwater  explosions  (3,  28-30).  To  predict 
the  effect  of  refraction  on  these  transients  at  caustics,  we 
have  combined  a  frequency  dependent  modified  ray  derivation 
that  accounts  for  refraction  (13)  with  a  Fourier  series  representa¬ 
tion  of  the  shock  wave  as  an  exponential  decay  from  a  peak 
pressure  (29).  With  this  approach,  we  have  successfully  matched 
experimental  pressure  records  measured  in  a  convergence  zone 
during  an  at-sea  test  involving  underwater  explosions  (8). 

However,  one  successful  comparison  does  not  validate  a  method. 

And  the  modified  ray  solution  used  to  account  for  refraction 
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effects  has  not  been  extensively  validated.  Therefore,  it  was 
felt  that  a  comparison  with  a  more  general  approach  was  necessary 
to  justify  our  continuing  use  of  modified  ray  results.  Normal 
mode  theory  seemed  a  logica.1  choice  for  our  standard  of  comparison 
since  it  is  more  general  than,  and  lacks  the  drawbacks  of, 
ray  and  modified  ray  theories. 

However,  combining  results  from  a  normal  mode  program  that 
accepts  an  arbitrary  sound  velocity  profile  with  a  transient 
source  representation  is  a  difficult,  time  consuming  task  that, 
to  this  author's  knowledge,  has  not  yet  been  successfully 
accomplished.  Furthermore,  if  ray  and  modified  ray  results  are 
adequate  for  underwater  explosion  shock  wave  predictions  near 
caustics  -  and  if  this  can  be  verified  without  complete 
comparisons  with  normal  mode  results  -  such  full  scale 
comparisons  may  be  unnecessary.  In  order  to  accomplish  this 
partial  validation,  we  decided  to  compare  modified  ray  theory 
with  normal  mode  theory  for  the  lower  frequencies  in  the  shock 
wave  (on  the  order  of  5°-100  Hz).  At  these  low  frequencies 
we  would  expect  discrepancies  between  normal  mode  theory  and 
modified  ray  theory  -  a  "high  frequency"  approximation  -  to  be 
most  apparent.  If,  on  the  other  hand,  they  agreed  fairly  well 
in  this  low  frequency  domain,  we  could  expect  them  to  agree  at 
least  as  well  at  the  higher  frequencies  in  the  shock  wave.  We 
could  then  have  confidence  in  the  use  of  modified  ray  theory 
on  all  parts  of  the  shock  wave  spectrum  without  resorting  to 
total  pulse  reconstruction  involving  normal  mode  theory.  We 
decided  to  do  these  comparisons  for  a  relatively  arbitrary  sound 
velocity  profile  so  that  we  would  be  testing  the  theories  for 
realistic  and  often  encountered  propagation  conditions. 
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These  comparisons  led  to  several  expected,  and  unexpected, 
chores.  The  finite  difference  normal  mode  program  (22)  required 
modification  in  order  that  it  work  properly  for  the  deep-water 
profiles  under  consideration.  A  mode  summing  program  had  to  be 
written.  The  ray  program  (6,  31)  being  used  had  to  be  made 
compatible  with  modified  ray  results  previously  derived  (13). 

And  a  ray  sorting  program  had  to  be  written  that  could  incorporate 
the  effects  of  one  or  more  caustics  into  the  intensity 
calculation. 

At  one  point  the  normal  mode  program  was  yielding  incorrect 
results  for  propagation  loss  versus  range.  We  found  ourselves 
using  ray  theory  and  modified  ray  theory  as  the  standards  of 
comparison,  with  normal  mode  theory  the  sample  requiring  validation 
a  complete  reversal  of  their  intended  roles.  We  will  discuss 
this  problem  in  order  to  point  out  some  of  the  difficulties 
associated  with  finite  difference  calculations  and  to  stress 
the  need  for  validation  of  any  model.  It  is  not  a  revelation 
to  state  that  any  propagation  model  (any  program,  for  that  matter) 
can  print  out  apparently  reasonable  propagation  loss  versus 
range  curves  while  doing  an  incorrect  calculation. 

In  the  following  sections  we  will  first  discuss  the 
common  roots  of  all  three  theories.  We  will  then  describe  the 
specific  ray,  modified  ray,  and  normal  mode  calculations  used 
in  these  comparisons.  After  checking  normal  mode  results  for 
a  shallow  water  model,  we  will  examine  a  deep  ocean  profile, 
concerning  ourselves  only  with  the  first  caustic  in  the 
convergence  zone  and  frequencies  of  50  Hz  and  100  Hz.  '.’hen  for 
a  slightli  shallower  profile,  we  will  consider  the  various 
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caustics  present  in  a  typical  convergence  zone.  Finally  we 
will  look  at  a  profile  yielding  several  more  nearly  horizontal 
caustics.  This  will  turn  out  to  be  the  most  difficult  to 
treat  with  the  horizontal  expansion  form  of  modified  ray 
theory  that  we  are  using,  and  here  we  will  discuss  a  useful 
alternative  derivation  (14). 

It  is  hoped  that  all  of  these  comparisons  will  give  a 
better  feel  for  the  strengths  and  weaknesses  of  the  various 
theories  as  well  as  indicate  where  each  one  can  be  used  to  best 
advantage.  They  also  should  help  explain  what  propagation 
paths  and  interference  mechanisms  are  causing  the  complex 
propagation  loss  curves  predicted  by  norma?,  mode  theory  in 
convergence  zones.  Finally  the  comparisons  will  demonstrate 
that  modified  ray  theory  is  valid  -  and  often  quite  useful  - 
in  predicting  pressures  on  and  near  caustics.  For  a  realistic 
deep  ocean  profile,  a  combination  of  ray  theory  and  modified  ray 
theory  results  often  yields  a  satisfactory  prediction  of 
pressures  in  the  convergence  zone  when  compared  to  a  more 
general  normal  mode  calculation.  The  fact  that  the  modified 
ray  theory  has  an  explicit  frequency  dependence  means  one  can 
locate  the  caustic  once,  calculate  a  few  parameters,  and  then 
find  the  intensity  for  any  frequency  rapidly,  This  is  as 
opposed  to  normal  mode  theory,  where  one  has  to  calculate 
intensity  at  a  given  point  independently  for  each  frequency. 

Thus,  modified  ray  theory  is  valid  for  the  treatment  of  the 
various  frequencies  in  the  shock  waves  from  underwater  explosions 
near  caustics,  as  well  as  being  an  inherently  easier,  more 
rapid  calculation  to  make. 
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2.  Wave  Equation  Roots  of  Ray  Theory,  Modified 
Ray  Theory  and  Normal  Mode  Theory 

As  we  have  said,  there  are  many  alternative  solutions  to 
the  wave  equation.  More  than  one  approach  leads  to  ray  theory 
type  results,  normal  mode  results,  or  any  other  category  of 
solutions.  Each  approach  emphasizes  a  particular  aspect  of 
propagation  phenomena.  In  this  section  we  will  summarize  the  two 
basic  derivations  that  have  served  as  the  basis  for  our  work 
in  modified  ray  theory  and  normal  mode  theory.  We  will  first 
go  through  the  details  of  the  Sachs  and  Silbiger  derivation  (13) 
that  arrives  at  a  ray  theory  solution,  chows  what  assumptions 
cause  it  to  break  down  at  caustics,  and  finally  arrives  at  a 
modified  ray  theory  solution  by  making  the  appropriate  changes 
in  the  derivation.  The  main  thrust  of  this  derivation  is  tc 
show  how  modified  ray  theory  springs  logically  from  ray  theory's 
failure  at  caustics  (ray  theory  itself  can  probably  be  better 
understood  by  examination  of  other  derivations  (1)).  Then  we 
will  start  again  with  the  wave  equation,  but  this  time  we  will 
go  through  the  details  of  Newman  and  I.ogenito's  (22)  finite 
difference  normal  mode  solution.  So  this  section  will  summarize 
the  basic  equations  from  previous  efforts  that  are  modified, 
evaluated,  and  compared  in  the  following  sections. 

We  start  with  the  reduced  wave  equation  for  pressure  and 
an  assumed  e  time  dependence! 

V*  +  kV\%( fc')  °-Po  ft*)  Tcfe)  (2.1) 

where  TUtf  =  CCtO/Cfi)  (2.2) 

1<  = 
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Using  cylindrical  symmetry,  Sachs  and  Silbiger  (13)  rewrite 
this  as» 


3*t>  +  ^4?  *fclnxa)p (2.4) 

,Sr'-  r  Or  7&r 


This  equation  is  then  multiplied  by  T6ttqr)  ,  and  integrated 


from  0  to  oo  to  obtaim 


+  K*  L  n*U)  -  S  -P«  fCt) 

att 


(2.5) 


(2.6) 


(2.?) 

(2.8) 


Where  *  $e  T6(K  $rWr  (2.6) 

and  Kr>*>=  (.  -ffS.llT.CKS^KVs  (2.7) 

•  iO«'»>l'"'tKSr)KVs  (2.8) 

Equation  (2.5)  is  then  solved  using  a  WKB  approximation  to 
arrive  ati 

a  [♦(*•>  ♦  ♦<*)] 

•fr  (c,2)  s  Jji_e _ (?  q) 

r  4V Mn'-SM''*  (‘-SMW  (2,9) 

.  ^  i/ 

where  )  r  \  3  d  %  ‘  i-.  is  the  turning  (2.10) 

J  point  depth. 

These  are  the  pertinent  equations  for  rays  beyond  their 

turning  points.  Equation  (2.8)  for  pressure  is  evaluated 

using  Equation  (2.9)  and  the  asymptotic  form  of  the  Hankel 

function  for  large  arguments i 


(2.9) 


p(r,  i)  a  Ps 
1  <HT 


n _ 1*  e.iktJ<!'r'il"^p 

Ci-S‘)'/lj  <2.n> 


W  t  J .  r,  *  I : 


(2.12) 
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To  use  a  stationary  phase  approach,  and  with  ray  theory  in 
mind,  Sachs  and  Silbiger  look  for  zeros  |t  of  in  the 

range  £  <  i,  hU.)  .  They  then  assume  that  the  main  contribu¬ 
tions  to  the  integral  are  in  the  vicinity  of  these  zeros. 

This  corresponds  to  ray  theory  in  that  (*•  will  turn  out  to 
be  Sin -fro  »  where  •G’q  is  the  initial  angle  of  a  ray,  As  long 
as  we  move  along  this  ray  path,  the  value  of  the  integral  will 
be  significant  and  energy  will  propagate. 

The  phase  term  W($)  is  then  expanded  around  these  zeros 

(13) » 

WC$)<  (2.13) 


where  W  {%<)  * 

and  w'cy-.o  by  definition  on  a  ray. 

Thus  Equation  (2.11)  becomes! 

PM=  5  f  if_ Hi _ _1  *. 

'  W  -*>  [L^TTr  C'-vVVlJ 

e*p  [-u 

Assuming  the  integral  in  Equation  (2,16)  is  slowly  varying 
when  £  is  near  and  dropping  higher  order  terms  than 

Vi"  t^)  yields  I  '/* 

P  (r,l)  -  £  | - -li - - - i  . 

exp  4 

This  is  a  ray  type  solution  for  the  refracted  wave  (beyond  the 
ray  turning  point),  where  each  term  in  the  sum  corresponds  to 
a  ray  leaving  the  source  at  angle  ('S’elj  1  Sin  with  the  vertical. 


(2.14) 

(2.15) 


(2.16) 


(2.17) 


k 
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We  note  that  the  amplitude  of  each  ray  is  determined  in  part 


=  [4"(2a)+4'U)] 


(2.18) 


where  Vl\)  was  the  only  derivative  of  the  Taylor  series 
expansion  (Equation  2.13)  assumed  to  contribute  to  the  integral 
(Equation  2.16). 

A  caustic  can  be  defined  in  several  ways!  One  alternative 
is  to  look  at  rays  at  a  constant  depth  and  plot  range,  r,  versus 
(Figure  2.1).  We  see  that  the  caustic  point  at  depth  z 
is  a  range  minimum  point  (it  could  also  be  a  range  maximum 
point).  Thus  for  the  caustic  point! 

i  =  sly  I  -.o 

'S*S«  (2.19) 

The  quantity  U\)  that  we  assumed  was  larger  than  all  higher 
derivatives  turns  out  to  be  zero  on  the  caustic.  Then 
Equation  (2.18;  becomes  infinite,  and  the  pressure  (Equation  2.17) 
becomes  infinite.  This  then  is  the  way  ray  theory  breaks  down 
at  the  caustic, 

It  is  a  natural  step  to  now  include  the  third  derivative 
in  the  integral  for  pressure.  However  Sac..s  and  Silbiger  go 
through  several  steps  in  inserting  the  derivative  that  are 


mathematically  correct,  but  not  at  all  obvious.  So  we  will 
explain  this  part  of  the  derivation  (32)  in  detail.  If  we  are 


on  a  caustic  at  (€6(l^  <  "the  ray  tnrough  the  caustic  point 

is  giv''”  by  •  Then 


1  ) 


(2.20) 


10 


NOLTH  lk-95 


This  is  just  the  criterion  for  being  on  a  ray  (Equation  2,15), 
applied  to  a  particular  caustic  ray.  Bj  definition  on  the 
caustic,  we  also  know 

W  (  8  0  (2.21) 


We  note  that  W  is  independent  of  range,  r,  since  r  vanishes 
from  Equation  (2.12)  upon  taking  two  derivatives  with  respect 
to  t.  „  Thus  if  we  are  at  depth  ,  but  at  a  range  r  not 

equal  to  rc,  we  can  be  on  a  ray,  but  not  at  the  caustic i 


,r,0  *°  W" 


(2.22) 

(2.23) 


However  if  we  replace  by  in  these  two  equations,  we  are 
not  on  a  ray  (the  ray  passes  through  rQ  at  zc),  sot 


But  since  there  is  a  caustic  at  depth  zc,  and  the  next 
derivative  is  independent  of  r  (r  or  r  ».akes  no  difference): 


W  t 

So  if  W  is  expanded  around  ^  for  points  off, 


(2.2 5; 


as 


well  as  on,  the  caustic,  then  W  will  be  zero  for  all  r, 


not  just  rc,  therefore  1 


,1  (2.26) 

We  note  that  is  non  zero  for  r  ^  rc.  Also  from 

Equation  (2.12): 

Tr+^**c«]-|V  (2.27) 
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Since  Ul$c  (the  caustic  ray  satisfies  the  ray 

criteria) 


VJ  (  $t, «  r- 


rt  s  AT 


(2.28) 


If  we  now  also  expand  the  amplitude  ^actors  in  Equation  (2.17) 
around  ^  ,  drop  higher  than  first  order  terms  as  leading 

to  higher  frequency  terms  after  integration,  and  put  in  the 
other  derivative  information  obtained  above,  we  get  Sachs  and 
Silbiger’s( 33) integral  expression  for  pressure  on  and  near  the 
caustic* 
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(2.30) 


and  (+)  in  s  is  used  depending  on  the  sign  of  w.  ,  The 
integral  is  xher.  expressed  in  terms  of  the  Airy  function 
(Figure  2.2),  yielding: 


h 


p(w 


*  T$  K 
4* 


LrcCnv.^)’/*0-sJ)T/t. 


Ai(t^)  «xp  OkWt-cn*)o 


(2.3D 


where  (  +  )  in  the  Airy  function  is  determined  by  (wfi  •<  °)  , 

Equation  (2.31)  is  valid  on  both  sides  of  the  caustic,  as  well 
as  on  the  caustic  at  Ar*0  ,  We  will  go  into  considerable 
detail  in  later  sections  as  to  this  equation’s  applicability 
under  various  circumstances.  In  general,  it  will  turn  out  to 


29) 
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be  quite  useful  in  predicting  pressures  on  the  caustic  itself 
as  well  as  in  the  adjacent  shadow  zone  and  caustic  region. 

And  the  Airy  function  will  be  important  in  understanding  how 
pressure  near  the  caustic  varies  with  frequency  and  distance 
from  the  caustic,  AT  . 

We  now  start  again  with  the  wave  equation  -  this  time  for 
velocity  potential  -  and  summarize  the  Newman  and  Ingenito 
derivation  (ZZ)  resulting  in  a  finite  difference  normal  mode 
solution.  They  start  off  with  the  Helmholtz  equation  in  both 
a  water  and  a  fluid  bottom  layer  (Figure  2.3) s 


and 


for  the  region  0  i  ?,  £  H 


5  0 

ct 

(2.33) 

for  the  region  H  J  i  i 

The  boundary  conditions  are 

(2.34) 

/“,  iM  = 

(2.35) 

^  U*H  U»H 

(2.36) 

the  conventional  pressure  release  surface,  and  continuity  of 
impedance  and  normal  particle  velocity  at  the  bottom  boundary. 
Equation  (2. 32)  leads  to  a  solution  for  ^  as  a  sum  of  discrete 
and  continuous  modes.  Discrete  modes  dominate  at  ranges 
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beyond  a  few  water  depths,  so  these  are  all  that  are  considered. 
(When  we  examine  the  node  sum  in  the  normal  mode  section,  it 
will  be  clear  why  contributions  from  modes  in  the  continuum  die 
out  with  range.)  Separating  variables,  and  normalizing  the 
depth  variable  ,  we  are  left  with  the  ^  dependent 

wave  equations  to  solve  (3*0* 


(2.3?) 


4*2*  W  ♦  H1, 

[j£.  -<] 
L  efa) 

<0 

2h  cn)  “>o 

0  i  ^  i  1 

dV^tt)  4  H1 

-  ,  « n 

5Lh  (*l)  •  o 

d*1 

1  £  *1  i  00 

(2.38) 


where  is  the  eigenvalue  for  the  nth  eigenfunction. 

Newman  and  Ingenito  divide  the  water  layer  into  m  equal  finite 
difference  layers,  each  with  velocity  C,^)  (Figure  2.4). 

This  is  accomplished  by  defining  m  layer  depths, y^  ( l/m  (...  ...|  , 

and  linearly  ir.terpo_  ting  between  the  input  sound  velocity 
points  to  find  e  velocity  at  each  of  these  points.  Then  into 
Equation  ( ? . 3V  '  n  the  water  layer,  they  substitute  a  first 
central  difference  for  the  second  derivative t 


_____  — (2.39) 

where  h  is  the  incremental  layer  depth  H/m 


Now  the  differential  equation  is  in  finite  difference  form* 
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Thus  we  now  have  an  expression  relating  the  wave  function 
amplitude  for  the  nin  mode  at  the  i+1  depth  point  to  the  amplitude 
at  the  i^*1  and  i-1  points,  the  velocity  at  the  i**1  point  and 
the  wave  number  kn<  We  will  discuss  in  Section  (5)  how  this 
is  used  to  find  the  individual  modes,  and  how  they  are  then 
summed  to  find  the  propagation  loss. 
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Figure  2.2  Airy  Function 
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Figure  2.3  Two  Fluid  Half-Space  Model 
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3.  Ray  Theory 

It  was  felt  that  rather  than  add  one  more  ray  tracing 
program  to  the  already  infinite  number  in  existence,  we  could 
modify  one  already  existing  to  suit  our  purposes.  We  wanted 
first  of  all  to  be  able  to  find  all  the  ray  paths  to  a  given 
point  in  space  and  add  the  arrivals  coherently  in  order  to 
calculate  the  net  intensity  at  that  point.  Then  we  wanted  to 
be  able  to  modify  the  program  to  calculate  the  actual  pressure 
on  the  caustic  using  modified  ray  theory  (Equation  (2.31)^. 
Finally,  we  wanted  to  be  able  to  add  in  shadow  zone  contributions 
at  a  point  of  interest,  when  that  point  is  in  the  shadow  zone 
of  one  or  more  caustics. 

We  chose  to  use  and  modify  CONGRATS  (Continuous  Gradient 
Ray  Tracing  System),  a  program  written  by  Weinberg  and  Cohen 
at  NUSC  (31).  CONGRATS  fits  the  velocity  profile  data  points 
with  layers  of  the  following  depth  dependence « 

c  -1  ',/i 

C(a)  r  )  V0+  + 

l  “Fft-wiV )  (3-1) 

By  the  appropriate  choices  of  the  four  arbitrary  parameters 
vQ,  gQ,  g-p  g2  (zQ  is  the  depth  at  the  top  of  the  layer), 
one  can  either  describe  the  profile  as  a  series  of  straight 
line  segments  each  with  the  appropriate  gradient  (Figure  3.1A), 
or  as  a  series  of  curves  defined  so  that  both  the  velocity  and 
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gradient  are  continuous  at  each  layer  interface  (Figure  3.1B).* 
Once  the  profile  is  specified,  travel  time,  range,  and  amplitude 
along  a  ray  can  be  easily  obtained  since  the  integrals  for 
these  quantities  are  readily  evaluated  in  terms  of  elementary 
functions  (31). 

CONGRATS  can  be  used  to  plot  ray  diagrams  in  order  to 
qualitatively  examine  the  sound  field  for  a  specific  profile 
(Figure  1.1).  But  it  can  also  be  used  to  quantitatively 
evaluate  the  intensity  at  specific  points  in  the  sound  field. 
This  is  accomplished  in  the  following  way.  One  specifies 
the  velocity  profile  and  defines  the  appropriate  profile  layers. 
Then  the  source  depth  is  defined,  as  are  the  particular  range- 
depth  target  points  of  interest.  Finally,  one  specifies  a 
grid  of  rays  (by  source  angle)  of  fine  enough  spacing  and 
sufficient  angular  width  in  order  to  ensonify  the  target  range- 
depth  points  with  each  significant  type  of  ray  arrival  (It 
is  no  smal]  task  to  pick  the  ray  spacing  and  grid,  and  skill 
at  this  ray  selection  increases  with  practice).  CONGRATS  then 
traces  one  ray  at  a  time  to  each  target  depth  of  interest. 

When  CONGRATS  finds  two  consecutive  rays  that  bracket  a 
target  range  at  the  target  depth  (Figure  3.2),  it  assumes 
that  a  ray  between  those  two  (called  an  eigenray)  would  reach 
the  target  point.  So  the  program  then  interpolates  between 

♦The  continuous  layers  can  eliminate  many  of  the  false  caustics 
caused  when  rays  turn  at  layer  interfaces  and  the  gradients  are 
discontinuous  (35).  In  the  cases  we  treated  this  was  not  a 
problem,  however,  and  constant  gradient  layers  were  adequate. 
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the  two  rays  in  order  to  find  the  travel  time,  amplitude,  etc., 
for  that  particular  eigenray.  It  then  was  written  to  add  all 
the  eigenrays  at  each  target  range-depth  point  coherently  or 
incoherently.  By  doing  a  series  of  range  points  at  a  fixed 
depth,  one  could  then  generate  a  typical  propagation  loss 
versus  range  curve  (Figure  3.3). 

As  originally  written,  CONGRATS'  sorting  programs  for 
adding  the  eigenrays  coherently  v/ere  not  compatible  with  the 
CDC  6400  computer.  Furthermore,  it  could  not  calculate 
actual  intensity  on  the  caustic  when  a  ray  passed  through  a 
caustic.  And  when  points  of  interest  were  in  the  shadow  zones 
of  one  or  more  caustics,  it  could  not  add  in  these  shadow 
zone  contributions  to  those  from  "real"  rays  passing  through 
these  points.  These  are  the  changes  made  at  the  Naval  Ordnance 
Laboratory  by  the  author  with  the  help  of  Jean  Goertner  and 
Robert  Thrun. 

In  this  oecuion,  we  will  discuss  the  coherent  sorting 
program  written  for  a  CDC  6400  computer.  This  program 
calculates  the  resultant  intensity  at  each  target  point  of 
interest.  It  also  can  add  in  shadow  zone  contributions  at 
each  point  for  any  number  of  caustics,  and  eliminate  bottom 
bounce  arrivals  if  this  is  desired.  In  the  modified  ray 
section,  we  will  discuss  how  CONGRATS  has  been  modified  to 
calculate  intensity  on  the  caustic  using  a  modified  ray 
calculation  (13). 

As  we  have  said,  CONGRATS  finds  the  eigenrays  that  pass 
through  each  target  point  by  interpolating  between  each  pair 
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of  rays  that  bracket  each  point.  If  we  specify  enough  rays  of 
different  type  (direct,  surface  reflected,  etc,),  to  ensonify 
the  region,  we  have  at  each  point  a  complete  set  of  eigenrays 
that  account  for  all  significant  energy  paths  between  source 
and  receiver.  CONGRATS  then  writes  on  tape  or  permanent  file 
several  blocks  of  data.  Each  case  -  defined  as  the  set  of 
input  data  cards  up  to  the  next  process  card  -  results  first 
in  the  writing  of  a  data  block  containing  all  of  the  input 
data  in  CONGRATS  internal  units  1  kyds,  kyds/sec,  seconds  and 
radian.  The  second,  and  if  necessary,  succeeding  blocks  of 
output  data  are  groups  of  100  arrivals,  each  an  eigen: ay.  At 
this  point  they  are  sorted  by  increasing  source  angle  with  range- 
depth  points  all  intermixed. 

Appendix  (I)  is  a  listing  of  the  program  written  at  NCL  to 
sort  these  arrivals  and  add  them  coherently.  Each  eigenray  of 
source  angle  G .  that  passes  through  a  target  point  (r.  z)  has  a 

vl 

travel  time  T.  (r,  z)  seconds,  propagation  loss  L.  (r,  z)  dB,  and 
J  J 

additional  phase  shift  (from  surface  reflections,  caustics,  etc.) 
of  |>,  (r,  z).  Ignoring  shadow  zones  and  caustics  for  the 
moment,  one  can  then  find  the  net  propagation  loss  at  a  given  point 
(r,  z)  and  a  specific  frequency  t*J  byi 

N^,;z)rao  f3o/(pt(r,i>|  ae  (3.2) 

is  Source  pressure  I 

(-L  (r,l)/aoJ  .  r  _  ,,  v“»\ 

~  pc  ^  lo  £xp(i  [u.T(rt?)4  4(rti)Jj  (3.3) 

Pj  (r,  z)  is  the  net  pressure  at  point  (r,z).  It  is 

obtained  by  the  coherent  sum  of  real  rays.  This  summation 

is  what  the  main  body  of  the  program  accomplishes.  In 

the  Drogram  listing  in  Appendix  (I),  lines  13-3'.'  read  the 
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input  data  block.  Lines  34-90  read  in  the  blocks  of  eigenrays 
into  a  doubly  subscripted  array  Buffer  (I,  J)  and  convert  from 
CONGRATS  units  (kyds,  kyds/sec,  etc.)  to  original  input  units 
(meters,  meters/sec,  etc.),  (feet,  feet/sec,  etc.)  or  any  other 
of  the  many  combinations  allowed  in  CONGRATS.  The  portion  of 
the  program  from  line  (133)  to  the  end  first  sorts  the  arrivals 
by  increasing  target  depth  and  then  by  increasing  target  range 
for  each  depth.  Then  all  the  arrivals  at  a  given  target  point 
(r,  z)  are  added  in  phase  according  to  Equation  (3.3). 

Appendix  (I),  pages  (  149-155)  show  first  the  output  data  block 
containing  all  the  original  data,  page  (  149  ),  then  the 

sorted  collection  of  arrivals,  pages  (150-153).  and  finally  the 
resultant  intensity  and  phase  at  each  target  depth-range  point, 
pages  (154-155). 

For  some  studies  it  is  necessary  to  eliminate  all  rays 

that  reflect  off  the  bottom  (or  surface).  For  example,  our 

normal  mode  calculations  have  been  done  (for  simplicity)  with 

a  bottom  fluid  whose  impedence  is  matched  to  the  v/ater  column 

( °2  =  c2  (z  =  H),^3 2  =  =  1).  So  we  would  not  expect  first 

ordpr  bottom  reflections  to  be  present  in  the  normal  mode 

calculation,  And  we  would  not  want  the  ray  calculations  to 

include  bottom  reflections  either.  The  easiest  way  to  eliminate 

bottom  reflections  is  to  not  include  rays  beyond  the  grazing 

ray  in  the  calculation.  But  sometimes  this  is  not  advisable,* 

*For  example,  CONGRATS  needs  at  least  one  ray  at  each  depth  range 
point  in  order  to  store  the  point  for  coherent  sorting.  Sometimes 
when  we  are  in  tne  shadow  zone,  there  are  no  real  direct  rays  - 
only  bottom  reflected  rays.  We  need  these  rays  to  "save"  the 
points  of  interest  so  that  we  can  add  in  the  shadow  zone  arrivals 
at  these  points  using  the  summing  program. 
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and  another  way  of  eliminating  bottom  reflec Lions  must  be  used. 
CONGRATS  has  the  capability  for  including  bottom  (or  surface) 
loss  as  a  function  of  incident  angle.  So  if  we  want  to  drop 
out  bottom  (or  surface)  reflections,  we  set  the  loss  for  all 
angles  of  reflection  to  be  9000  dB.  Then  in  the  sorting  program, 
Appendix  (I),  page  (14?),  lines  (173-I84),  ve  chuck  to  see  if 
the  propagation  loss  for  any  arrival  is  greater  than  9000  dB 
(i.e.  it  has  reflected  off  the  bottom  at  least  once).  All  of 
these  arrivals  are  then  eliminated  from  the  calculation. 

As  we  have  said,  the  sorting  program  must  also  be  able  to 
add  in  shadow  zone  arrivals  whenever  the  target  point  is  in 
the  shadow  zone  of  a  caustic.  V/e  have  chosen  to  work  with  a 
modified  ray  expression  that  yields  results  at  various  ranges 
on  either  side  of  the  caustic  (as  well  as  on  it)  at  the  fixed 
depth  the  caustic  occurs  at.  A  ray  of  source  angle  Qq  may  have 
a  caustic  at  range  rc  and  depth  zQ .  Then  the  contribution 
from  this  caustic  at  any  range  r  at  the  caustic  depth  z£  is 
given  by  1  i/2 


pCr' LnMV(,...J  *;tv>  o-»> 

(t:  vj;  0) 

where  all  quantities  of  interest  were  defined  in  Section  (2), 

Equa  ;"ns  (2.30-2.31).  The  pressure  expression,  and  its 
region  of  validity,  are  further  discussed  in  the  modified  ray 
section.  For  negative  p  and  1^1  >  1  we  are  in  the  double  arrival 
region  associated  with  the  caustic.  In  Section  (4),  we  will 
discuss  the  use  of  the  caustic  solution  in  this  region  as 
opposed  to  the  actual  two  arrivals  as  calculated  from  ray  theory. 
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However,  in  the  shadow  zone  there  is  no  problem,  Ray 
theory  yields  no  shadow  zone  contribution,  and  the  modified 
ray  theory  contribution  is  necessary  for  completeness.  We 
need  to  evaluate  both  the  amplitude  and  the  phase  of  these 
arrivals,  so  that  we  may  incorporate  them  into  the  coherent 
sum  (Equation  (3*2)).  When  a  ray  passes  through  a  caustic, 
CONGRATS  prints  out  the  amplitude  on  the  caustic.  This  is 
Equation  (3.^)  with  Ar  =  0,  so^=  0.  This  is  essentially! 

pirc^t)  s  (3 .3) 

since  the  remainder  of  the  amplitude  expression  is  the  same  on 
and  off  the  caustic.  Off  the  caustic  at  arbitrary  r,  the  pressure 
is  then 

r  B  A!  re'll-  BA'.C  K^Ar]  (3.6) 


So  the  pressure  off  the  caustic,  in  terms  of  the  pressure  on 
the  caustic  calculated  by  CONGRATS  isi 


|>CrM2c)  lA;(f)l/A;£o) 


(3.7) 


And  the  propagation  loss  at  distance  Ur  off  the  caustic  is 

=  *ao  Lo^(rCl2t)  I k\(p)\  j Ai(o)j  (; 

or 

aoU^Jj>(.r  *t)]  :  =  AtfPDB  -  |A.'(p)|  /. ( 

where  AMPDB  is  the  propagation  loss  on  the  caustic  calculated 
by  CONGRATS.  We  also  need  the  total  phase  of  the  arrival. 


We  note  that 

Wt(r)  =  §tr  *  +  <j>(^ 


(3.10) 
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It  is  defined  in  such  a  way  (Equation  (2.10))  that  4(z0) 
and  4(z_)  must  be  integrated  analytically  or  numerically. 

C 

But  we  would  like  to  use  quantities  already  calculated  by 

CONGRATS  where  possible  to  avoid  unnecessary  calculations. 

We  can  follow  in  reverse  a  derivation  by  Officer  (36)  in 

order  to  arrive  at  a  relationship  between  W.(r  )  and  T(r,)t 

the  travel  time  to  the  caustic  alrec.dy  calculated  in  CONGRATS. 

We  take  W ( r  ) , 

c  c  _ 


w ,(rt)  5  $  r„ 

'A  '/, 

*  {  Cn*U)-^]da.  [rt'c*)-§l3  d* 

(3.11) 

*  J 

1T 

C*4  "A 

*  ^  da 

(3.12) 

*0 

where  £  -  Sin(^4)  t  CCft^SlfVd* 

CCa) 

(3.13) 

,  vk 

•-  •*c(a4)\  [**(*)  -C^)x]  da 

(3.1*0 

2* 

where  ^  ~  S.in(i8^) 

C.U)  C(a) 

(3.15) 

Therefore: 

Wc.(rc) 

C(%.) 

s  ^  [**(*>  -  C$V]  *"cU 

(3.16) 

But  for  a  ray: 

3Vlt  (v* ) 

i  0  z  Y'  *  ^  4(2o)  +  "S  4(2  ") 

(3.17) 
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Them 


r  s  -  'S  -  a<j>ta) 

»s 


(3.18) 


or  taking  the  derivative  of  Equation  (3.12), 


r  =  i  ^  /l  da.  (3.19) 


Substituting  this  into  Equation  (3.16)  for  r=rc  and  combining 
terms i 


CClo) 


J  L  L^V).(s;)1]4 


da  (3.20) 


[  «\i) 

l  L  [«*«-($•]* 


da 


(3.21) 


*t 


CU) 


(3.22) 


but  1'  =  iiD_±t  =>>  £0S^t  [l- (§')lCaU)] '/i(3.23) 


Therefore  Vk^X) 

CUo) 


da 


C(4)  CoS-b-t 


ii.  --  T  (Vt)  (3<24) 


cu> 
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So  if  we  have  travel  time  to  the  caustic  from  CONGRATS,  T(rc), 


we  can  immediately  evaluate  W(rc)  byi 

S  c£%»)  T(rc)  (3.25) 

Or  for  arbitrary  r,  off  the  caustic  in  either  direction, 

VittM  S  {jr  +  4^*3*  4^‘^  (3.26) 

s  *  4f*#)  +  4(s<^ +  icr* (3-27) 

s  |tCr.rt>  •♦VJcCr,.)  (3.28) 

*-  |tAr  4CCa.)T(rc)  (3.29) 

Then  HWt(r)  1  ^  St  T*OV)  (3.30) 

CUO 


Thus  the  main  phase  term  in  Equation  (3.14)  is  readily  calculated 

from  constants  ( ^ c  ,  W  ,  C ( ) ,  6r)  and  the  travel  time  to 

the  caustic  point  T(rc)  calculated  by  CONGRATS.  To  complete 

the  phase  we  must  also  add  “tr/4  and  any  extra  phase  shift  <|(n) 

associated  with  extra  caustics,  surface  reflections,  etc. 

the  ray  may  have  undergone.  We  can  then  add  to  the  arrivals 

given  by  Equation  (3.3)  shadow  zone  arrivals  for  the  caustics 

# 

occurring  at  the  depth  of  interest  (sc)  and  various  ranges 
rQ(l),  rc(2)  ...  rc(n)« 


(*%[>  .Vl/ao) 


P,p  is  the  coherent  sum  of  all  rays,  including  caustic  shadow  zone 
contributions.  This  is  done  in  the  sorting  program,  Appendix  (I), 

4*  h 

lines  (93-130).  We  only  add  an  arrival  from  *-’ne  r.  caustic  as 

long  as  we  are  in  the  shadow  zone  of  that  particular  caustic  (±fl>o). 

s  ' 

For  a  typical  range  minimum  ccX/erg  zone  caustic,  W^" 

29  X 

\ 

X 


HOLTS  7*»-95 


is  negative,  so  we  use  Ai(- p  ).  Thus  A  r  must  be  negative 
in  order  for  to  be  positive,  and  r  must  be  less  than  rc. 

So  we  are  to  the  left  of  the  caustic,  nearer  the  source. 

This  then  describes  CONGRATS  and  the  sorting  program  we 
have  written  to  accompany  it.  The  sorting  program  adds  the 
regular  ray  arrivals  coherently.  It  adds  in  shadow  zone  arrivrls 
coherently  when  the  target  point  is  in  the  shadow  zone  of 
one  or  more  caustics.  And  bottom  reflections  can  be  eliminated 
when  this  is  desired. 
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4.  Modified  Ray  Theory 

While  ray  theory  is  usea  extensively  for  propagation  loss 
calculation,  its  failure  near  caustics  is  well  documented.  The 
predictions  of  infinite  peak  pressure  on  the  caustic  and  zero 
pressure  in  the  adjacent  shadow  zone  are  two  unrealistic  and 
unacceptable  conditions.  This  failure  occurs  because  of  the 
"high  frequency"  approximation  inherent  in  the  eikonal 
equation  (37).  It  yields  a  picture  in  which  all  energy  within 
a  ray  tube  of  a  given  cross  section  at  the  source  remains  within 
that  ray  tube  as  the  cross  section  goes  to  zero  at  the  caustic 
(defined  as  the  locus  of  points  where  infinitesimally  close 
rays  cross).  So  the  finite  amount  of  energy,  contained  in  a 
zero  cross  section  ray  tube,  yields  infinite  peak  pressure. 

If  we  are  interested  in  pressure  on  the  caustic,  or  in  the 
adjacent  shadow  zone,  another  solution  must  be  used. 

To  this  end,  several  authors  (12-14)  have  solved  the 
wave  equation  for  expressions  that  remain  valid  on  the  caustic. 
Sachs  and  Silbiger  (13)  ana  Brekhovskikh  (12)  obtained  an 
expression  that  yields  pressure  on  the  caustic,  as  well  as  off 
the  caustic  horizontally  in  either  direction  (Figure  4,1a). 
Ludwig's  derivation  (14)  yields  pressure  on  the  caustic,  and 
off  the  caustic  in  either  direction  normal  to  it  (Figure  4. IB). 
Ludwig's  solution  has  both  a  uniform  and  non-uniform  result. 

The  uniform  asymptotic  theory  is  exact  in  the  ray  double  arrival 
region,  while  the  non-uniform  solution  is  only  approximately 
correct  there.  Ludwig's  non-uniform  result  is  very  similar  to 
Sachs  and  Silbiger's  result  (13).  The  similarities  and 
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differences  between  the  various  solutions  will  be  discussed  in 
more  detail  later  on  in  this  section.  The  horizontal  expansions 
in  particular  have  been  used  to  predict  pressure  histories 
from  underwater  explosions  in  convergence  zones  (8)  and  in  flooded 
quarries  (8,9).  They  also  have  been  compared  with  ray  theorv 
and  normal  mode  theory  for  a  monotomically  decreasing  sound 
velocity  profile  (6). 

In  our  earlier  work  (8),  we  used  a  constant  gradient 
profile  for  which  the  modified  ray  solution  could  be  obtained 
analytically.  We  will  start  out  by  summarizing  this  derivation 
for  pressure  on  and  near  a  caustic.  Then  we  will  discuss  the 
evaluation  of  the  various  quantities  necessary  to  actually 
calculate  pressure.  And  finally  we  will  relate  these  expressions 
to  the  appropriate  quantities  in  CONGRATS  that  were  being,  or 
are  now  being,  calculated. 

As  we  have  pointed  out  in  Section  (2),  one  can  solve  the 
wave  equation  using  a  WKB  approximation  and  arrive  at  a  ray 
solution  (Section  2,  Equation  2.1?).  This  equation  breaks  down 
at  a  caustic.  A  further  derivation  results  in  a  new  integral 
(2.29)  for  pressure  on  and  near  a  caustic.  When  the  integral 
is  evaluated,  Sachs  and  Silbiger  (13)  arrive  at  an  expression 
valid  at  caustics  (2.31)i 

Vi 

«■ 

KC  *t)  =  Po  - a."  (i/»>  exp  (t'KVJe-ijt  )  (4.1) 

<TT  [r,  (.n%.  * 


where  =  sin(0Q)t,  ( 0^ )c  i s  the  source  angle  passing  through 
a  caustic  at  (rc,  zc)  and  the  remaining  quantities  are  defined 
in  Equation  (2,30). 
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For  a  typical  convergence  zone  caustic  (Figure  4.2),  Wcm(rt) 

is  negative}  hence  we  use  (-fi).  Then  for  r  >  rc,  to  the 

right  of  the  caustic  at  f  =  0,  the  Airy  function,  and  so 

the  pressure,  first  rises  as  we  move  through  the  caustic 

region  (Figure  2.2).  Further  to  the  right  (  Vi  > '  ) ,  we  see 

oscillations  in  the  Airy  function  that  are  typical  of  the 

well  known  ray  double  arrival  region  adjacent  to  the  caustic. 

To  the  left  of  the  caustic,  for  positive  arguments  of  the  Airy 

function,  the  pressure  decays  exponentially  with  distance  off 

the  caustic  as  well  as  frequency  to  the  two  thirds  power. 

f!l 

For  a  caustic  with  Wc  {rt)  positive,  the  picture  is 
completely  reversed  with  the  shadow  zone  to  the  right  of  the 
caustic,  etc. 

Several  quantities  must  be  evaluated  in  order  to  calculate 
pressure  on  or  near  a  caustic  using  Equation  (4,1).  First  we 
have  to  know  when  a  particular  ray  goes  through  a  caustic. 

As  has  been  pointed  out.  a  caustic  is  a  range  minimum  or  maximum 
point  at  constant  depth.  So  at  a  <.austic 


£r 


-  £1*^  / j>JL  \  s  o 
WCv  / 


(4,2) 


where  Sm*Oj  s  £  s  C(30} 

0  v 

and  Cy  is  the  sound  velocity  at  the  vertex  of  the  ray. 


along  a  ray,  and 


CONGRATS  determines  this  by  evaluating 

3CV 

checking  for  cnanges  in  sign  of  this  quantity,  indicating  it 


has  passed  through  zero.  This  is  equivalent  to  the  check  on 


we  have  used  in  previous  Work  (38),  So  for  a  given 


source  angle  (0  )  we  can  move  along  the  ray  until 

0  TJCv 
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approaches  zero.  At  this  point  (r„,  z„)  we  are  on  the  caustic. 

By  tracing  a  number  of  rays  of  different  source  angles,  we  can 

generate  one  or  more  caustic  curves  in  r,  z  space  and  find  the 

caustics  at  a  particular  depth  of  interest  (Figure  4.3).  Then 

r  ,  z  ,  Ar  (to  the  point  of  interest),  n  (z  ),  C  ,  k,  and 
c  c  c  J  c 

Wc(r)  (Section  2,  Equations  2.10  and  2.12)  are  fixed  for  each 
caustic. 

Thus  we  need  only  to  evaluate  )>  (Equation  2.30)  in 
order  to  find  the  pressure  on  the  caustic.  We  have  an  expression 
for  £  in  Sachs  and  Silbiger  notation,  but  we  first  need  the 
proper  derivation  in  CONGRATS  notation.  In  CONGRATS,  the 
range  increment  in  a  layer  is  defined  ast 


irl  a  AR 


(4.3) 


Because  C(z)  is  defined  analytically  as  a  function  of  z 
(Equation  3*1) »  this  integral  can  be  evaluated  in  terms  of 
elementary  functions  (39).  Its  first  derivative,  'dr/<)Cv  , 
which  is  needed  for  ray  amplitude  calculation  and  caustic  location, 
was  also  evaluated  analytically.  But  we  need  the  second 
derivative ,  ’aV/'DCv  ,  for  caustic  amplitude  calculation. 

This  derivative  was  obtained,  and  it  is  summarized  in 


Appendix  (IV).  It  was  then  inserted  into  CONGRATS.  In  order 


to  use  this  derivative,  we  must  relate  it  to 


»WcW 


the  derivative  in  our  modified  ray  notation  (Equation  2.30). 
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We  start  off  withi 


Wfr(r^  “  §  r  +  4(*0)  ♦  4^*^ 
*^VJfc(r)  »  f  +  ^  4Ut)  4-  *3  4(4«) 


(4.4) 

(4.5) 


On  a  ray 


From  Equation(4. 5) . 


But  from  Equation(4. 


Then 


-^4(40 
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(4.6) 

T>i  1 

ssst 

^HJ  , 

7>  §*• 

s  i  r  . 

3$  L  -j§ 

. 
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./^vJ 

as  long 
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V 

as  we  are  on  e 

i  ray. 

i 

;  .ir 

(4.8) 

where  r  is  the  analytic  range  expression  in  CONGRATS. 


Since 

'i)*'  -  <"2f" 

*  ■Sc«*  H 

(4.9) 

-  -'W  s  r  T>r  ^cj\ 

(4,10) 

-  -1  { ^ 

«>$  ^  ^  '  acv  'H1, 

(4.11) 

But 


chf  r  0 


on  a  caustic,  so 
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So  with  the  second  range  derivative  (Appendix  IV)  we  have 

added  to  the  program,  and  the  appropriate  constants, 

3  ^ 

we  can  calculate  d  Vldrt)j2f  and  so  &  .  This  then  is  the 
final  quantity  necessary  to  calculate  pressure  on  the  caustic 
by  Equation  (4.1),  and  evaluate  modified  ray  theory  arrivals  in 
Equation  (3*31)* 

Now  that  we  have  sho vm  how  to  calculate  the  pressure  on 
and  near  a  caustic,  we  would  like  to  discuss  the  validity  of 
its  use  at  various  locations  near  the  caustic  and  its 
relationship  to  other  methods  of  caustic  calculation.  Looking 
at  Equation  (4,1),  we  can  picture  the  pressure  as  a  constant 
times  the  Airy  function  (Figure  2.2).  For  Ai(0),  we  are  on 
the  caustic  and  the  expression  is  valid  for  most  single,  well- 
behaved  caustics.  Fcr  Ai(~p  ),  the  pressure  at  first  grows 
in  what  we  call  the  caustic  region  (0  £.  \p  \  £\ ,S  ).  Here 
we  expect  the  caustic  solution  to  be  better  than  the  ray 
solution  which  is  diverging  (Figure  4,4).  We  note  that 
f>  =  k2^t  fcr.  The  quantity  t  is  roughly  tne  focusing 
factor,  indicating  the  strength  of  the  caustic.  It  is  also 
related  to  the  slope  of  the  caustic.  For  typical  convergence 
zone  caustics  we  characterize  as  vertical  (Figure  4.2)  (depth 
to  range  slope t  lilO),  6  is  on  the  order  of  .002.  Therefore 
the  source  frequency  and  width  of  the  caustic  region  (where 
ray  theory  is  invalid)  are  related  byt 


Ar  £  150 /\^'l 


(4,14) 
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Thus  a  100  Hz  source  v/ill  have  a  caustic  region  about  2500 

meters  wide.  The  lower  the  frequency,  the  wider  tne  caustic 

region  where  ray  theory  is  not  valid. 

For  larger  negative  arguments  (to  the  right  or  left  of 

in.  . 

the  caustic  depending  on  whether  W  (rt)  is  negative  or 

positive),  we  are  in  the  ray  double  arrival  region,  where 
ray  theory  is  valid.  The  caustic  solution  (with  an  asymptotic 
form  of  the  Airy  function)  does  yield  two  arrivals  (40): 


However,  they  are  of  equal  amplitude,  a  condition  that  we  only 
expect  to  be  true  very  near  the  caustic. 

Since  this  expression  is  only  an  approximation  derived 
from  the  value  on  the  caustic,  we  shall  see  that  from  case  to 
case  the  oscillation  pattern  agrees  more  or  less  with  the 
actual  pressure  as  calculated  by  ray  theory  or  normal  mode 
theory.  This  agreement  and  disagreement  can  be  understood  by 
considering  Figure  (4.5).  Consider  a  caustic  at  point  (A) 
(Figure  4.5A).  We  then  extend  the  caustic  solution  horizontally 
into  the  double  arrival  region  to  point  (B).  This  extension 
will  only  be  valid  when  the  rays  actually  passing  through 
point  (B)  have  essentially  the  same  history  near  the  caustic 
as  the  ray  passing  through  the  caustic  at  point  (A).  Only 
this  way  can  the  ray  going  through  (A)  "know"  what  rays  at  (B) 
should  look  like. 

In 
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For  the  simple  circular  caustic,  this  is  easiest  to  see 
(Figure  4.5A).  As  long  as  the  radius  of  curvature  is  roughly 
the  same,  rays  at  (B)  (having  passed  through  C)  are  roughly 
the  same  as  the  ray  passing  through  a  caustic  at  (A).  Then 
an  extrapolation  to  (B)  based  on  (A)  is  roughly  correct.  This 
is  easiest  to  see  at  Point  (E)  when  ray  (1)  is  at  a  distance 
from  the  caustic  equal  to  AB,  The  ray  pattern  looks  similar 
to  that  at  (B).  But  if  the  curvature  changes  (Figure  4.5B) 
then  the  rays  at  (B)  (having  passed  through  C)  do  not  have  the 
same  history  as  the  ray  passing  through  (A).  Again  look  at 
(E),  when  ray  (1)  is  the  rppropriate  distance  from  the  caustic. 
The  ray  pattern  at  (E)  is  different  than  at  (B),  while  the 
expansion  assumes  the  pattern  is  the  same  at  both  points. 

For  more  complicated  caustics,  it  is  not  as  easy  to  see.  But 
the  same  rule  of  thumb  applies.  As  long  as  the  radius  of 
curvature,  or  slope,  cf  the  caustic  does  not  change  much  along 
the  caustic,  say  from  (  c)  to  (A),  we  can  extend  the  caustic 
solution  horizontally  to  a  point  near  (A)  through  which  rays 
from  (C)  pass. 

In  the  shadow  zone,  we  have  no  real  rays  from  that  caustic 
to  compare  modified  ray  theory  to.  The  shadow  zone  arrival 
is  the  only  one.  But  from  the  Airy  function  (Figure  2,2), 
the  pressure  falls  as  the  two  thirds  power  of  frequency  at  a 
given  distance  hr,  or  linearly  with  horizontal  distance  from 
the  caustic  for  a  fixed  frequency.  We  can  only  verify  this 
behavior  by  comparison  with  normal  mode  theory,  since  ray 
theory  predicts  zero  pressure  in  the  shadow  zone. 
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For  more  nearly  horizontal  caustics  (Figure  4.6)  that 
have  depth  to  range  slopes  of  approximately  It  100,  ^  is  less 
than  ,001.  In  this  case  the  width  of  the  caustic  region  is 
roughly  two  to  five  times  wider  than  in  the  vertical  case 
(Equation  4.16).  Intuitively  we  would  expect  the  peak 
amplitude  to  be  lower,  as  the  energy  is  spread  over  a  broader 
caustic  region.  This  is  generally  the  case.  The  horizontal 
cau-tics  tend  to  have  significantly  lower  peak  amplitude, 
with  the  pressure  falling  off  much  more  slowly  into  the  shadow 
zone  than  in  the  vertical  case. 

The  same  arguments  hold  for  extensions  into  the  double 
arrival  region  from  smooth  horizontal  caustics  that  hold  for 
vertical  caustics.  Reasonable  predictions  depend  on  the 
caustic  having  the  same  radius  of  curvature  throughout  the 
region  of  interest.  However  for  horizontal  caustics  (Figure  4.6) 
the  rays  must  travel  considerably  farther  to  reach  an  arbitrary 
point  (B);  distance  AB  from  the  caustic.  So  the  caustic  must 
maintain  the  same  local  slope,  or  radius  of  curvature,  for  a 
much  longer  distance.  This  makes  it  more  difficult  for  the 
horizontal  expansion  to  work,  and  it  is  this  case  where  a 
normal  expansion  is  more  reasonable.  In  Figure  (6),  FB  is 
much  shorter  than  AB.  Thus  a  prediction  for  (B)  based  on  (F) 
is  more  likely  to  work  than  one  from  (A)  extended  to  (B). 
Furthermore,  there  may  not  be  a  caustic  at  the  proper  depth, 
line  (D),  In  this  case,  a  horizontal  expansion  into  the  shadow 
zone  is  not  possible. 
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Figure  4.4  Airy  Function  and  Divergent  Ray  Solution 
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5.  Normal  Mode  Theory  and  Program  Validation 

Normal  mode  theory  is  the  most  general  of  the  three 
solutions  to  the  wave  equation  that  we  have  examined.  By 
solving  the  wave  equation  directly,  rather  than  the  ray 
equation  approximation,  we  avoid  the  problems  inherent  in 
ray  theory  and  modified  ray  theory.  Normal  mode  theory  does 
not  break  down  at  caustics.  It  yields  finite  non-zero  pressure 
in  the  shadow  zone,  accounting  for  energy  diffracted  from  the 
caustic.  And  it  yields  the  resultant  intensity  at  a  given  point 
directly.  No  ray  theory  type  addition  of  various  ray  energy 
paths  -  and  the  worry  of  missing  one  -  need  be  done. 

This  is  not  to  say  that  normal  mode  theory  is  flawless. 

In  solving  for  the  net  intensity,  it  eliminates  one  of  che 
nice  things  about  ray  theory  -  the  very  physical  picture  the 
ray  diagram  (Figure  4.2)  gives  about  the  paths  that  energy 
is  taking  between  source  and  receiver.  There  have  been  attempts 
(41)  to  relate  various  portions  of  the  mode  spectrum  to  specific 
types  of  rays.  However,  we  are  more  comfortable  understanding 
the  intensity  in  terms  of  the  ray  $.nd  modified  ray  arrivals 
interfering  with  each  other,  and  then  comparing  this  to  the 
normal  mode  theory  result  in  order  to  understand  the  significance 
of  the  various  intensity  variations  in  the  propagation  loss 
curve.  Other  faults  of  normal  mode  theory  will  be  discussed 
during  the  remainder  of  this  section.  But  in  general,  for 
deep  ocean  profiles  it  becomes  a  case  of  solving  for  many  modes 
(which  at  some  point  becomes  too  many),  and  doing  it  completely 
over  each  time  a  different  source  frequency  is  specified.  The 
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frequency  dependence  is  not  like  ray  or  modified  ray  theory, 
where  geometry  is  considered  first  as  ray  paths  are  calculated, 
a  id  then  frequency  is  accounted  for  only  during  the  coherent 
addition  of  arrivals.  In  normal  mode  theory  each  frequency 
results  in  a  completely  different  set  of  modes,  and  so  a 
completely  new  calculation  must  be  done  from  start  to  finish. 

So  as  with  anything  else,  each  theory  has  its  appropriate 
use.  Normal  mode  theory  is  suitable  for  shallow  water  (few 
modes) ^profiles  and  for  some  deep  water  (moderate  number  of 
modes)  profiles.  It  serves  as  an  excellent  standard  of  comparison 
with  each  to  check  out  ray  and  modified  ray  theory  results. 

Ray  theory  is  especially  useful  for  many  deep  water  profiles, 
and  modified  ray  theory  extends  its  usefulness  by  adding  the 
capability  for  intensity  calculations  on  and  near  caustics. 

Once  or.„  decides  that  a  normal  mode  solution  is  desired, 
there  are  several  alternative  routes  for  doing  the  calculation. 

As  discussed  in  Section  (2),  we  have  chosen  to  use  a  finite 
difference  approach,  enabling  us  to  deal  with  an  arbitrary 
velocity  profile.  We  have  taken  a  program  originally  written 
for  the  shallow  water  case,  and  we  have  made  the  appropriate 
modifications  so  that  it  could  handle  the  deep  water  case.  In 
the  process  we  have  learned  quite  a  bit  not  only  about  normal 
mode  calculations  in  general,  but  also  about  the  peculiarities 
of  finite  difference  calculations  in  particular. 

First  we  will  briefly  describe  the  use  of  the  equations 
obtained  by  Newman  and  Ingenito  (22)  to  calculate  the  normal 
modes  for  a  given  profile.  Then  we  will  discuss  the  summing 
program  written  to  calculate  the  propagation  loss  once  the 
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normal  modes  have  been  obtained.  Finally  we  will  discuss  the 
changes  made  in  the  normal  mode  computer  program  so  that  it 
could  treat  the  deep  water  case, 

As  discussed  in  Section  (2),  Newman  and  Ingenito  arrived 
at  a  finite  difference  expression  for  the  unnormalized  wave 
function  at  the  i  +  1  point  ?n  terms  of  the  sound  ,relocity  in 
the  i*"  layer,  the  trial  mode  wave  number  kn,  and  the  value  of 
the  v/ave  function  at  the  i  -  1  and  i  ^  ^  pointi 


The  velocity  profile  has  been  split  up  into  m  equal  depth  layers 
of  depth  h,  and  the  velocity  (c,).  specified  in  each  layer. 

If  we  had  two  points  to  start  off  Equation  (5.1).  Z 

n ,  i 

and  Z^P  ,  and  a  k  ,  we  could  generate  the  mode  shape.  Then 
n  i  4  n 

if  the  mode  shape  satisfied  all  the  boundary  conditions  and 
had  the  correct  shape  (the  n'  mode  has  n-1  zero  crossings), 
we  would  know  that  we  had  the  correct  value  of  kn,  as  well 
as  the  proper  mode  shape.  Using  the  two  boundary  conditions 
at  the  water-bottom  interface,  Newman  and  Ingenito  (42)  arrive 
at  expressions  for  the  first  two  wave  function  values,  at  the 
water-bottom  interface  and  one  layer  upi 


'  n.i 


-  Nfi 

■for  or  ^;| 

(5.2) 

Inh  AVi 

-  iPd 

-h’w  -kV1d  ••••  ] 

(5.3) 

l  * 

1 

X  j 

1 

■for  ^  :  l  -  V> 
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where 


o*. 

w  - 


H1 

Hl 


1 


(5.4) 


(There  are  several  typographical  errors  in  this  part  of  their 
report  and  this  part  of  the  program  listing  as  well.) 


The  program  then  takes  these  two  points,  the  trial  kn, 

Equation  (5.1).  and  the  velocity  profile;  it  follows  the 
wave  function  up  to  the  surface.  kn  is  iterated  on  until  we 
have  the  proper  number  of  zero  crossings,  and  the  wave  function 
is  within  some  small  epsiion  of  zero  at  the  surface.  This 
last  requirement  satisfies  the  last  boundary  condition  of  the 
surface  as  a  pressure  release  boundary. 

The  iteration  on  kn  should  be  started  in  a  way  that 
minimizes  the  number  of  trial  solutions  tested  before  the 
proper  kn  is  obtained.  For  the  first  mode,  a  maximum  mode 
wave  number  kn(max),  and  a  minimum  mode  wave  number  kn  (min)  are 
defined  byi 


kn  (max) 
kn  (min) 


u; 

LC't  >]  min 

W 


(5.5) 

(5.6) 


where  [cw  )]  m^n  is  the  minimum  sound  velocity  in  the  water 
and  C,  is  the  velocity  in  the  bottom.  Several  trial  kn 
between  kn  (max)  and  kn  (min)  are  tried  until  a  mode  shape  with 
no  zero  crossing  is  obtained*.  Then  this  kn  and  kn  (max) 
are  used  as  the  bracketing  values  which  the  final  iteration 


*The  first  mode  just  reaches  zero  amplitude  at  the  surface. 
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x.  u. 

process  starts  off  with.  For  the  n  mode,  kn_-^  is  used  as  a 
lower  limit  instead  of  kn(max),  and  the  upper  limit  is  adjusted 
down  from  kn(min)  until  n-1  zero  crossings  are  obtained.  Other¬ 
wise  the  iteration  orocess  is  the  same. 

The  normalization  constant  is  found  by  numerical  integration 
of  the  wave  function  in  the  water  layer  anu  analytical  integration 
of  it  in  the  bottom  layer.  After  each  normalized  mode  U  A*l) 
is  calculated,  the  wave  number,  k  ,  and  other  necessary  information 
about  the  mode  is  written  on  a  permanent  file  when  our  revised, 
version  is  used.  Also  written  on  this  file  is  the  particular 
mode  amplitude  at  one  source  depth  Un(  )  and  up  to  five 
receiver  depths  U  )  (Main  program,  lines  190-230).  The  file 
then  serves  as  input  for  the  summing  program  listed  in 
Appendix  (III). 

In  th*  s  summing  program,  v/e  calculate  the  pressure  at  a  given 
point  in  space  by  relating  pressure  to  velocity  potential  1 

P  ^  ’  (5.7) 

—  iidt 

Using  the  suppressed  e  time  dependence  and  the  velocity 
potential  (43) 1 


We  obtain  from  (5-7) 

*4'  (.Tr)'k  [  ‘•U’WU" 

The  continuous  portion  of  the  mode  spectrum  consists  of  modes  for 
which  J  knJ  <  4)/c2.  Solutions  of  Equations  (2.37)  and  (2.38) 
lead  to  imaginary  values  of  and  so  an  exp(-knr)  in 
Equation  (5-9).  So  these  modes  are  dam.ped  out  for  ranges  greater 
than  a  few  water  depths. 
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We  then  assume  spherical  spreading  close  in,  so  that  the 
reference  pressure  at  unit  distance  is  (44)i 

t  k  R 

5  £_  (5.10) 

■MTT  R 

which  at  1  meter  yields t 

tO£>|  nro  (5.11) 

Therefore  the  propagation  loss  isi 


Thus  given  the  mode  wave  numbers,  other  constants  such  as 
density  and  water  depth,  H,  and  mode  amplitudes  at  the  source 
and  receiver  depths  for  each  mode,  we  can  calculate  the  propagation 
loss.  This  is  then  what  the  summing  program,  Appendix  (III), 
does.  Since  there  is  nothing  special  in  Equation  (5.12) 
about  source  and  receiver  depth,  we  have  written  the  program 
so  that  any  depth  specified  in  the  normal  mode  program  can  be 
used  as  the  source  depth  in  tne  summing  program.  The  program 
reads  in  information  from  the  normal  mode  program  (Appendix  II), 
checks  what  source  and  receiver  depths  are  required,  calculates 
propagation  loss  according  to  Equation  (5.12),  and  plots  or 
prints  propagation  loss  versus  range. 

When  Pedersen  and  Gordon  (6)  compared  normal  mode  theory 
to  ray  theory,  they  pointed  out  that  this  interchangeability 
of  source  and  receiver  is  not  true  with  ray  theory,  only  normal 
mode  theory.  They  then  proceeded  to  show  that  the  ray  theory 
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result  should  he  modified  by  the  ratio  of  source  to  receiver 
sound  velocity  to  make  it  equivalent  to  and  comparable  to 
normal  mode  theory* 

Hnm  =  Hr  +  10  Lo®  (Cs/CK}  (5.13) 

While  the  correction  makes  ray  theory  more  "exact,"  it  is 
not  significant  for  many  realistic  profiles.  For  example, 
in  our  profile,  Figure  (4.2),  the  worst  case  would  be  source 
at  the  velocity  minimum,  receiver  on  the  bottom.  The  correction 
would  still  only  be  .1  dB.  The  correction  is  important  when 
some  analytical  profiles  are  used.  Here  the  receiver  or  source 
sound  velocity  may  get  very  low  (admittedly  unrealistic),  and 
the  correction  becomes  significant. 

Once  the  normal  mode  summing  program  was  written  and  checked 
out,  we  did  a  few  straightforward  comparisons  in  order  to 
verify  the  output  of  the  normal  mode  program.  We  compared 
results  from  our  program  to  results  from  two  other  programs 
for  a  shallow  water  profile  (Figure  5.1).  DiNapeli  (24) 
had  originally  used  this  profile  to  compare  his  fast  field 
program  (FFP)  to  Bartberger's  normal  mode  program  (23).  The 
FFP  uses  a  completely  different  approach  by  first  fitting  the 
sound  velocity  profile  with  exponential  layers.  It  then  solves 
the  wave  equation  directly  in  terms  of  Greens  functions  and 
uses  a  fast  fourier  technique  on  the  computer. 

Figure  (5.2)  shows  the  two  original  results  in  addition 
to  a  calculation  for  the  same  profile  using  our  nox'mal  mode 
program.  Discrepancies  in  the  original  comparison  were 
attributed  to  FFP's  use  of  exponential  layers  as  opposed  to 
Bartberger's  linear  layer  fit  of  the  velocity  profile  (Figure  5.1). 
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We  would  expect  our  results  to  be  closer  to  the  other  normal 
mode  result,  since  both  use  constant  gradient  segments. 

However  our  normal  mode  calculation  includes  only  "trapped 
modes"  (i.e.  modes  whose  phase  velocity  is  less  than  the  sound 
velocity  in  the  bottom  fluid).  This  is  the  discrete  portion 
of  the  mode  spectrum.  Bartberger's  normal  mode  program 
continues  calculating  modes  up  into  the  continuum  (phase  velocity 
greater  than  bottom  fluid  sound  velocity),  and  so  it  should 
yield  a  different  result.  But  the  comparisons  in  Figure  (5.2) 
show  essentially  the  same  propagation  loss  versus  range  for 
all  three  theories.  This  comparison  gave  us  confidence  in 
the  program  as  originally  written  ior  shallow  water  cases. 

We  then  proceeded  to  test  the  normal  mode  program  for  a 
deep  water  case.  We  chose  a  typical  deep  water  case  yielding 
a  well  developed  convergence  zone  (Figure  5.3).  This  ray 
diagram  for  a  305  m  source  depth  indicates  that  to  the  right 
of  the  caustic  (Figure  5«3.  line  AB)  bordering  the  convergence 
zone,  there  is  a  well  developed  double  arrival  region.  One 
arrival  has  passed  through  rhe  caustic,  the  other  is  approaching 
it.  We  then  ran  CONGRATS  for  a  500  m  receiver  and  100  Hz 
source  frequency.  We  chose  ranges  from  5^  km  to  6 2  km,  which 
took  the  receiver  from  the  shadow  zone,  over  the  caustic,  and 
into  the  double  arrival  region.  In  the  shadow  zone,  CONGRATS 
found  no  real  rays  as  expected  (bottom  bounce  arrivals  were 
dropped  from  the  calculation  as  discussed  in  Section  3).  To 
the  right  of  the  caustic,  two  arrivals  were  found  at  each 
range  point.  They  were  then  added  coherently  using  the  ray 
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summing  program  (Appendix  I).  A  Tty 2  phase  shift  was  added 
in  CONGRATS  to  the  arrival  that  had  passed  through  the  caustic. 
This  phase  shift  is  a  well  known  but  sometimes  controversial 
feature  of  ray  theory.  Figure  (5.4)  shows  the  propagation  loss* 
versus  range  curve  for  a  500  m  receiver.  The  propagation  loss  is 
essentially  infinite  to  the  left  of  the  caustic  -  which  is 
indicated  by  the  vertical  line  at  56.3  km.  The  propagation  loss 
jumps  to  zero  at  the  caustic  (infinite  intensity),  then 
increases  to  a  plateau  at  about  57  km.  From  this  point  on  the 
propagation  loss  curve  oscillates  as  the  two  arrivals  interfere 
first  destructively,  then  constructively.  Both  are  very  roughly 
equal  in  amplitude.  That  plus  the  Tt/Z  phase  shift  results  in 
This  particular  curve  (Figure  5.4).  There  have  been  many 
supportive  papers  (45-4?)  on  the  use  of  the  TT/2  phase  shift 
at  caustics,  and  later  we  will  add  our  own  example  of  how  this 
is  the  phase  shift  necessary  for  ray  theory  to  be  consistent 
with  modified  ray  theory  and  normal  mode  theory. 

In  any  event,  now  that  we  had  the  ray  theory  result,  we 
ran  the  normal  mode  program  for  the  same  profile.  We  split  the 
profile  into  the  maximum  allowable  number  of  layers,  1000,  and 
found  from  a  preliminary  run  that  there  were  139  modes.  To 
match  the  ray  calculation  geometry  we  set  the  bottom  density 
equal  to  1  and  the  bottom  fluid  velocity  equal  to  the  velocity 
in  the  wa  o  -t  the  bottom  (this  is  a  matched  impedance  with  no 
*A11  propagation  loss  values  will  be  in  dBtre  1  yd  unless 
otherwise  specified. 
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first  order  bottom  reflections).  Y/e  then  used  the  normal  mode 
program  to  calculate  all  the  mode  shapes,  stored  the  mode 
amplitudes  for  each  mode  for  source  (305  m)  and  receiver 
(500  m)  depths,  and  found  the  propagation  loss  using  the  mode 
summing  program  (Appendix  II;').  Figure  (5.5)  shows  this  curve 
as  well  as  the  ray  theory  curve,  They  couldn't  be  more  out  of 
phase  if  we  had  tried  to  make  them  so. 

This  disagreement  was  unexpected,  to  say  the  least.  We 
proceeded  to  examine  both  programs  in  order  to  find  the  problem 
and  evexvtuaily  found  the  difficulty  to  be  with  the  normal  mode 
calculation.  We  were  calculating  139  modes  and  using  1000 
finite  difference  layers  in  the  water  column.  This  meant  that 
for  the  higher  modes  that  were  oscillating  one  hundred  times 
or  more  in  the  water  column,  we  we re  allowing  10  finite  difference 
’’ayers  or  less  to  fit  each  oscillation.  The  scarcity  of  layers 
would  result  in  a  poor  representation  of  higher  modes  as  we 
followed  the  wave  function  from  one  layer  to  the  next,  and  more 
layers  would  enable  us  to  more  accurately  find  that  mode  shape. 

This  inaccuracy  in  mode  shape  would  affect  not  only  the  normalization 
constant  but  also  the  choice  of  mode  wave  number,  kn> 

Figure  (5.6)  shows  the  normal  mode  calculation  for  2000 
layers.  The  comparison  is  better,  but  still  not  satisfying. 

Figure  (5.7)  shows  the  normal  mode  calculation  for  3900  layers. 

Here  there  is  still  closer  agreement  with  ray  theory.  And  the 
change  in  the  normal  mode  calculation  from  2000  to  3000  layers 
is  small  compared  to  the  change  from  1000  to  2000  layers.  It 
was  felt  that  calculations  for  more  than  J000  layers  were 
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unnecessary,  and  a  rule  of  thumb  was  adopted  that  at  least 
20  times  the  number  of  modes  equaled  the  number  of  finite 
difference  layers  required  for  the  calculation. 

The  larger  number  of  layers  required  necessitated  a 
rewriting  of  the  program.  We  now  calculate  one  mode  at  a 
time,  extract  the  required  information,  and  go  on  to  the  next. 

This  is  as  opposed  to  the  original  storage  of  a  group  of  twelve 
modes  at  a  time.  The  program  was  substantially  rewritten  to 
keep  its  size  small  while  increasing  the  number  of  finite 
difference  layers.  A  listing  is  shown  in  Appendix  (II). 

We  also  made  several  changes  to  improve  the  accuracy  and 
speed  of  the  calculation,  using  knowledge  gained  while  looking 
for  the  fla w  in  the  program.  For  example,  regula  falsi  (48) 
is  usually  faster  than  simple  halving  as  a  root  finding  technique, 
and  as  such  it  was  used  exclusively  in  the  original  program. 

But  for  the  lower  modes  in  the  deep  v/ater  case,  we  found  that 
the  first  trial  values  of  kR  were  yielding  values  of  the 
unnormalized  wave  function  at  the  surface  on  the  order  of 
+  10100.  This  was  instead  of  the  very  small  value,  £<  1, 
needed  to  satisfy  the  pressure  release  boundary  condition. 

Because  of  this,  regula  falsi  was  taking  excessively  long  to 
zero  in  on  the  proper  kn.  So  we  now  use  halving  of  the 
difference  between  successive  kn's  until  the  wave  function  at 
the  surface  is  less  than  10"’°  (Appendix  II,  Subroutine  Half, 
lines  45-70),  and  then  we  use  regula  falsi  until  the  wave 
function  at  the  surface  is  less  than  6  .  In  general,  this 

cut  the  number  of  iterations  required  by  more  than  half  for 
the  lower  modes. 
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Another  change  in  the  program  involved  the  method  used 
for  zeroing  the  wave  function  near  the  surface.  Ideally, 
for  each  mode  one  can  find  a  kR  so  that  the  mode  has  the 
proper  number  of  zero  crossings  (n-1)  and  reaches  an  amplitude 
value  of  zero  at  the  surface.  In  practice,  using  numerical 
evaluation  of  the  mode  shape,  the  amplitude  is  never  quite 
zero  at  the  surface.  Even  for  changes  in  kn  of  one  pare  in 
10  ,  the  surface  wave  function  value  may  still  be  small,  but 

finite.  So  the  program  was  set  up  to  satisfy  the  surface 
boundary  condition  as  follows  1 

zn  (^  =  o)  i  0  (5.14) 

where  £  is  some  small  number  on  the  order  of  .01  or  less. 

If  kn  is  changing  by  less  than  one  part  in  lO^0,  and  the 
amplitude  still  isn't  within  €  of  zero  at  the  surface,  the 
iteration  process  stops.  At  this  point  the  wave  function  will 
look  like  Figure  (5.8A)  or  Figure  (5.8B),  It  will  approach 
zero  amplitude  at  the  surface,  miss  reaching  zero  amplitude 
at  the  surface,  or  cross  too  early  resulting  in  one  too  many 
zero  crossings.  This  information  hear  t/iS  surface  is  incorrect, 
and  essentially  useless.  Newman  and  Ingenito  (22)  decided  to 
zero  the  mode  amplitude  at  these  depths  in  the  following  way. 

They  calculated  the  mode  shape,  again  starting  from  the  bottom. 
They  checked  until  there  were  n-1  zero  crossings  (for  the  n^*1 
mode),  Point  A,  Figure  (5.8a).  Then  they  followed  the  wave 
function  as  it  first  increases  (AE)  and  then  decreases  beyond  B. 
At  this  point,  beyond  B,  the  wave  function  should  monotonically 
decrease  to  zero  at  the  surface.  If  it  does  (Figure  5.8A), 
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it  is  left  alone.  But  if  it  doesn't,  at  some  point  it  may  cross 
the  axis  or  start  to  increase  again  (Figure  5.8B).  This  is  the 
point  where  the  wave  function  is  assumed  to  misbehave,  and  it 
is  zeroed  from  there  to  the  surface.  This  is  done  in  subroutine 
(Iterate)  when  Flag  equals  2. 

While  this  approach  works  well  for  many  profiles,  it  can 
yield  incorrect  results  for  some  modes  with  some  profiles. 

For  example,  consider  an  asymmetric,  double  well  profile 
(Figure  5.9).  While  most  modes  present  no  problem,  consider 
mode  4  (Figure  5-10).  This  is  the  unnormalized  wave  function. 
The  program  would  follow  it  to  point  (A),  continue  up  to  (B), 
and  then  down  to  (C).  There  would  be  the  proper  number  of 
zero  crossings,  so  it  wouldn't  expect  the  amplitude  to  increase 
beyond  point  (B),  It  would  then  incorrectly  zero  out  the  mode 
beyond  (C)  in  the  upper  well  and  normalize  incorrectly, 
resulting  in  an  incorrect  mode  (Figure  5.11).  To  avoid  this, 
we  chose  to  zero  the  amplitude  starting  from  the  surface. 

The  possibilities  are  the  same  as  before  (Figure  5.8)  with 
the  first  case  being  left  alone,  while  the  other  two  require 
zeroing.  So  we  start  from  the  top  and  see  which  shape  the 
first  oscillation  has.  This  is,  of  course,  after  first 
generating  and  storing  the  entire  mode  shape  by  starting  at 
the  bottom-water  interface,  and  using  the  last  value  of  k 

n 

obtained.  This  mode  shape  is  not  correct  near  the  surface, 
but  it  is  the  best  we  can  do  within  the  accuracy  limitations 
of  the  method  and  computer.  We  then  zero  the  mode  amplitude 


62 


NOLTR  7*1-95 


near  the  surface  where  necessary.  Figure  (5.12)  shows  mode  4 
again,  hut  this  time  after  proper  zeroing,  w'  see  the  correct 
normalized  mode  shape.  These  changes  were  made  in  Subroutine 
Iterate,  Appendix  (II). 
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5.1  Environmental  Description  for  a  Single 
Exponential  Layer(From  DiNapoli( 24) ) 
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Figure  5-6  Profile  I,  Normal  Mode  (2000  inc.)  and  Ray  Theory  (z=500m,  f=100Hz) 
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Figure  5-7  Profile  I,  Normal  Mode  (3000  inc.)  and  Ray  Theory  (z=500m,  f=100Hz) 
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Figure  5*12  Mode  4,  Correctly  Normalized 

(•Small  Amplitude  Oscillation  in  Bottom 
Well  Roughly  One  Tenth  as  Large  as 
Shown  Here) 
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6,  Profile  I,  Comparisons  Near  a  Single  Caustic 

Our  first  comparisons  of  the  three  theories  were  done  using 
the  deep  ocean  profile  used  in  our  checkout  of  the  normal  mode 
program  (Section  5.  pages  57-60).  Profile  I*(Figure  5»3) 
was  specifically  chosen  so  that  the  caustic  bordering  the 
convergence  zone  (Figure  5.3.  line  AB)  would  have  a  well 
developed  double  arrival  region  to  the  right  of  it.  While  we 
are  also  interested  in  comparisons  close  in  and  for  more  than 
one  caustic,  we  felt  it  would  first  be  advisable  to  compare  the 
theories  in  the  simplest  region  where  all  three  are  valid, 
near  a  single  caustic.  To  the  left  of  the  caustic  (line  AB), 
there  is  at  first  a  simple  caustic-related  shadow  zone.  This 
region  will  not  turn  out  to  be  completely  free  from  effects  of 
the  surface  reflected  shadow  zone  boundary  close  in  (Figure  5.3* 
line  CD).  But  the  two  boundaries  are  far  enough  apart  so  that 
we  can  separate  out  the  effects  of  each.  To  the  right  of  the 
caustic,  because  of  the  bottom  depth  of  the  profile,  we  have 
a  double  arrival  region  free  from  other  arrivals  for  an 
appreciable  distance  from  the  caustic.  This  enables  us  to 
compare  the  three  theories  not  only  in  the  caustic  region,  but 
for  several  oscillations  in  the  double  arrival  region  as  well. 

Calculations  were  done  for  a  source  depth  of  305  m,  and 

receiver  depths  of  2$0  m,  500  m,  and  1500  m.  This  way  we  were 

treating  different  situations  -  a  receiver  shallower  than  the 

source,  between  the  source  and  sound  channel  axis  depths,  and 

near  the  sound  channel  axis,  Figure  (6.1)  shows  comparisons 

between  ray  theory  and  modified  ray  theory  at  the  two  deeper 
*Tnput  sound  velocity  profile  data  in  Appendix  V. 
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depths  for  a  source  frequency  of  100  Hz.  As  expected,  ray 
theory  propagation  loss  diverges  as  the  caustic  is  approached 
from  the  right,  and  reaches  zero  (infinite  peak  pressure)  at 
the  caustic  -  indicated  by  the  vertical  line  in  each  part  of 
the  figure.  Modified  ray  theory  indicates  a  finite,  non-zero 
propagation  loss  at  the  caustic,  and  the  intensity  falls  off 
exponentially  from  there  into  the  shadow  zone.  Just  to  the 
right  of  the  caustic,  for  Airy  function  argument  |yO  I  £  1,5, 
we  are  in  the  caustic  region  -  or  caustic  boundary  layer  (13). 
Here  we  expect  modified  ray  theory  to  be  better  than  ray  theory. 
This  is  the  region  in  each  part  of  Figure  (6.1)  where  the 
modified  ray  theory  propagation  loss  first  decreases  and  finally 
reaches  a  minimum,  before  again  starting  to  increase.  We  see 
that  at  about  this  minimum  point,  ray  and  modified  ray  theories 
merge  and  reach  their  best  agreement  (Figure  (6.1A)  at  5?  km, 
for  example). 

Beyond  this  point  we  expect  ray  theory  to  be  valid j 
modified  ray  theory  may  agree  with  it,  but  this  will  depend 
on  the  local  geometry  of  the  caustic.  At  500  m  (Figure  6.1A), 
modified  ray  theory  does  do  well  in  this  double  arrival  region 
and  follows  the  ray  theory  oscillations  adequately.  At  1500  m 
(Figure  6. IB),  however,  modified  ray  results  fail  rapidly  in 
the  double  arrival  region.  It  is  soon  almost  completely  out  of 
phase  with  the  ray  theory  results.  This  disagreement  is  not 
unexpected  (49),  and  points  out  that  care  must  be  exercised  in 
extending  modified  ray  theory  far  into  the  double  arrival  region. 
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In  Figure  (6.2),  we  now  add  normal  mode  theory  results  at 
each  depth.  On  the  caustic  at  each  depth,  modified  ray  results 
agree  with  the  normal  mode  results.  Thus,  modified  ray  theory 
does  yield  valid  predictions  on  the  caustic.  Then  looking  at 
500  m,  we  see  that  normal  mode  values  duplicate  the  modified 
ray  results  in  the  caustic  region,  and  then  it  agrees  with  both 
ray  aad  modified  ray  theories  in  the  double  arrival  region.  To 
the  left  of  the  caustic,  the  normal  mode  results  agree  with  modi¬ 
fied  ray  results  in  predicting  an  exponential  drop  off  in 
intensity  into  the  shadow  zone.  Only  when  the  pressure  has 
fallen  35  to  40  dB  and  a  weaker  diffraction  pattern  from  the 
close  in  boundary  and  bottom  starts  to  predominate,  does 
modified  ray  theory  fail.  So  in  this  case,  it  predicts  the 
caustic  related  shadow  zone  effect  as  far  as  we  can  see  it. 

Now  we  look  at  1500  m  where  ray  and  modified  ray  disagree 
in  the  double  arrival  region  (Figure  6. IB).  Since  we  are  using 
normal  mode  theory  as  the  standard  of  comparison  valid  everywhere, 
we  expect  it  to  agree  with  modified  ray  theory  in  the  caustic 
region  (0  £  i  -1.5)  -  where  modified  ray  is  valid.  Then  we 

4 

expect  normal  mode  results  to  agree  with  ray  theory  results 
in  the  double  arrival  region  (beyond  52.5  km)  -  where  ray 
theory  is  valid.  Figure  (6.2B)  verifies  this  behavior.  In 
the  shadow  zone,  modified  ray  results  are  still  in  good  agreement 
with  normal  mode  results.  Thus  the  tendency  of  modified  ray  to 
be  valid  in  some  region  to  the  right  of  the  caustic  appears 
to  be  accompanied  by  a  similar  tendency  to  the  left  of  the 
caustic.  Sachs  (49)  has  explored  this  problem  of  the  shadow 
zone  validity  of  modified  ray  theory  using  an  idealized  model 
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resulting  in  a  circular  caustic.  He  demonstrated  that  as  one 
proceeds  into  the  shadow  zone  (large  positive  arguments  p  of 
the  Airy  function),  modified  ray  theory  breaks  down  when  compared 
to  an  exact  solution  of  his  particular  problem  (Figure  6.3).  So 
if  we  could  follow  modified  ray  farther  into  the  shadow  zone  - 
and  we  will  do  this  later  in  a  multi-caustic  case-  we  would  expect 
it  to  fail.  However,  for  our  realistic  profile,  modified  ray 
theory  is  still  in  good  agreement  with  normal  mode  theory  when 
the  intensity  is  quite  low. 

Sachs  (49)  also  considers  the  use  of  complex  ray  theory 
in  the  shadow  zone.  Complex  ray  theory  is  the  ray  type  solution 
valid  in  shadow  zones  where  regular  ray  theory  is  not  useful. 

From  Figure  (6.3),  we  can  see  that  complex  ray  theory  does  work 
well  in  the  shadow  zone.  It  is  somewhat  more  difficult  to  use, 
since  it  requires  solutions  for  complex  roots  of  the  ray  equations 
for  range,  time,  etc.  However,  it  can  be  a  useful  addition  to 
ray  solutions  for  work  in  the  deep  shadow  zone. 

We  next  did  the  same  comparison  for  50  Hz  (Figure  6.4).* 

The  agreement  among  the  three  theories  is  about  the  same  as 
for  ICO  Hz.  There  are  fewer  oscillations  over  the  same  range 
increment  because  of  the  lower  frequency.  Furthermore,  the 
disagreement  between  ray  and  modified  ray  theory  at  1500  m 
(Figure  6,4c)  starts  to  occur  at  a  different  range  -  at  about 
56  to  57  km  -  between  the  first  and  second  nulls.  V/hat  is 
roughly  constant  at  the  point  of  disagreement  is  p  ,  the 
argument  of  the  Airy  function  -  which  is  a  function  of  the 
frequency  through  k2//^  and  also  distance  off  the  caustic,  h  r. 


*Here  we  include  a  250  m  receiver  for  the  first  time. 
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Thus  for  the  same  geometry  caustics  (constant  a  ),  and  different 
frequencies,  f  is  the  factor  that  scales  the  Airy  function  by 
spreading  or  shrinking  it  over  the  space  coordinate  r.  Finally, 
Figure  (6,5)  shows  the  results  for  50  and  100  Hz  in  the  shadow 
zone  at  500  m  depth.  Modified  ray  predicts  that  intensity  will 
fall  off  exponentially  v/ith  and  distance  off  the  caustic, 

Ar,  Figure  (6.5)  demonstrates  that  both  of  these  modified  ray 
theory  approximations  are  good  estimates  of  the  behavior 
predicted  near  the  caustic  by  normal  mode  theory.  Finally,  it 
should  also  be  noted  that  the  consistency  of  all  three  theories 
in  the  double  arrival  region  depends  on  the  ir/2  phase  shift  that 
is  inserted  ir.  the  ray  theory  result  for  the  ray  that  has  passed 
through  the  caustic.  This  phase  shift  is  also  implicit  in  the 
modified  ray  theory  result  (50).  Only  this  way  do  they  agree 
with  normal  mode  theory.  So  this  is  just  one  more  demonstration 
of  the  presence  of  a  'ir/2  phase  shift  ?t  a  caustic. 
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Figure  6.1A  Profile  I,  Ray  and  Modified  Ray  Theory 
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Fieure  6.2A  Profile  I,  Normal  Mode,  Ray,  and  Modified  Ray  Theories 

Cs  =500m,  f=100Hz) 
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Figure  6,4a  Profile  I  Normal  Mode,  Ray,  and  Modified  Ray  Theory  (z=250m,  f=50Hz) 
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Figure  6.5  Frequency  Dependence  of  Normal  Mode  and  Modified 
Ray  Theory  in  Shadow  Zone  Adjacent  to  the  Causti 
(z=500m,  f=50,  100Hz) 
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7.  Profile  II,  Comparisons  for  a  Multi-Caustic  Convergence  Zone 

While  our  main  concern  is  with  caustics  and  the  convergence 
zone,  for  completeness  we  compared  ray  theory  to  normal  mode 
theory  in  the  near  field  region  to  the  left  of  the  shadow 
zone.  Figure  (7.1)  shows  the  velocity  prof ile* and  ray  diagram 
considered.  Figure  (7.2)  shows  the  near  field  comparisons  for 
the  three  depths  of  interest.  Because  we  are  using  the  asymptotic 
expression  for  the  Hankel  function,  good  for  large  arguments,  we 
do  not  expect  the  normal  mode  result  to  be  valid  within  a 
range  of  about  one  water  depth  of  the  source.  And  this 
breakdown  does  show  up  (for  example,  at  ranges  less  than 
10  km  in  Figure  (7.2C))  where  the  normal  mode  result  starts  to 
diverge. 

In  Figure  (7.2),  the  propagation  loss  curve  oscillates 
several  times  as  the  direct  and  surface  reflected  rays  interfere. 
Then  at  a  specific  range  at  each  depth  (29.6  km  at  15u0  m, 
for  example),  the  ray  theory  propagation  loss  increases  to 
infinity  as  we  enter  the  shadow  zone.  Normal  mode  results 
continue  to  indicate  a  finite  amount  of  energy  present  in  the 
shadow  zone.  At  250  m  (Figure  7.2A),  the  general  fall  off  in 
diffracted  energy  with  increasing  range  is  most  obvious. 
Superimposed  on  this  fall  off  is  an  oscillation  due  to  energy 
that  is  probably  reflected  by  the  matched  impedance,  discontinuous 
gradient  bottom.  As  the  receiver  depth  becomes  deeper,  the 
distance  between  the  near  shadow  zone  boundary  and  caustic 
shadow  zone  boundary  decreases  and  the  bottom  grows  closer, 
so  that  by  the  tims  we  reach  1500  m,  the  average  propagation  loss 
remains  roughly  constant  throughout  the  shadow,  zone. 

*Input  sound  velocity  profile  data  in  Appendix  V. 

89 


IIOLTR  7h~95 


We  next  examined  the  convergence  zone  resulting  from  the 
same  profile  (Figure  7.1).  The  various  caustics  present  in 
the  ray  diagram  are  shown  in  Figure  (7.3).  The  receiver 
depths  rre  the  same  as  before,  250  m,  5°0  m,  and  1500  m. 

However  due  to  the  complexity  of  the  ray  pattern,  we  will 
discuss  each  receiver  depth  separately. 

Table  (7.1)  summarizes  the  rays  passing  through  each 
range  point  at  a  depth  of  250  m  in  the  convergence  zone. 

Figure  (7.'+)  shows  the  ray  theory  propagation  loss  calculated 

N 

by  CONCRATS  which  summed  coherently  the  rays  present  at  each 
point,  it  also  contains  the  propagation  loss  versus  range 
cur\'e  as  calculated  by  normal  mode  theory.  In  the  ray  theory 
curve,  we  see  a  caustic  at  58.9  km.  Two  arrivals  interfere 
to  the  right,  of  this  up  to  point  (  U  ) ,  where  the  bottom  cuts 
off  one  of  them.  The  single  arrival  region  extends  to  point  (V), 
where  we  pass  over  the  next  caustic  and  pick  up  one  more  arrival. 
This  caustic  is  branch  B  in  Figure  (7.3) •  It  is  the  surface 
reflected  branch  made  up-  of  those  single  arrivals  adjacent  to 
the  first  caustic  that  were  not  cut  off  by  the  bottom.  At 
point  (W)  we  pass  over  another  caustic  resulting  from  rays 
that  are  reflected  off  the  surface  near  the  source  (Figure  7.3, 
Branch  C).  Again  we  pick  up  another  ray.  Finally  at  point  (X). 
and  beyond,  we  pick  up  another  arrival  that  has  reflected  off 
the  surface,  passed  through  a  caustic  (Figure  7.3,  Branch  H),  and 
reflected  off  the  surface  again. 

There  are  two  causes  for  abrupt,  discontinuous  changes  in 
the  propagation  loss  curve  calculated  by  ray  theory.  First, 
the  appearance  of  the  caustic  usually  results  in  a  sharp  spike 
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in  the  propagation  loss  curve.  This  is  seen  at  58.9  km,  at 
the  border  of  the  convergence  zone.  The  modified  ray  results 
we  add  in  will  smooth  these  out  as  shown  in  the  previous 
comparisons  (Figure  6,1),  The  other  cause  of  abrupt  changes 
is  the  sudden  cut  off  of  a  given  group  of  rays  by  the  bottom. 

This  problem,  also  discussed  by  Pedersen  and  Gordon  (16),  is 
present  in  our  matched  impedance  bottom.  In  any  case,  a 
particular  series  of  rays,  each  of  which  turns  deeper  and  reaches 
a  different  range  point  in  a  smooth  pattern,  suddenly  is  cut 
off  by  the  bottom  (Figure  ?.5).  There  is  not  a  smooth 
transition  as  rays  beyond  the  critical  ray  are  cut  off  and  so 
do  not  continue  the  previously  established  smooth  pattern. 

Norma],  mode  results  do  indicate  a  smooth  transition,  and  this 
problem  of  ray  theory  we  have  not  treated. 

We  do,  however,  treat  the  caustic  problem.  We  add  in 
modified  ray  results  whenever  v/e  are  in  the  shadow  zone  or 
on  the  caustic.  Furthermore,  we  add  in  the  modified  ray 
results,  and  take  out  the  two  divergent  ray  arrivals,  whenever 
we  are  in  the  caustic  region,  0>p>-1.5in  the  Airy 
function. ^  Figure  (7.6)  shows  the 'combination  of  ray  theory 
and  modified  ray  theory,  as  compared  to  normal  mode  theory. 

As  in  other  comparisons,  the  convergence  zone  boundary  is 
adequately  described  by  modified  ray  results.  At  point  (Y), 
we  see  the  start  of  an  oscillation  caused  by  the  single  "real" 
ray  interfering  with  a  shadow  zone  arrival  from  the  caustic 

*Wherev  r  we  talk  of  putting  in  modified  ray  results  in  the 
caustic  region,  it  will  be  understood  that  we  take  out  of 
the  sum  of  all  rays  those  particular  rays  associated 
with  that  caustic. 
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at  point  (  V).  In  general,  ray  plus  modified  ray  goes  a  long 
way  toward  matching  the  normal  mode  results.  As  discussed 
previously,  the  worst  points  of  comparison  are  where  the 
bottom  cuts  off  an  arrival,  Points  (U)  and  (Z).  Normal  mode 
results  indicate  a  smooth  transition,  but  simple  ray  theory 
has  no  way  of  treating  this.  Instead  an  arrival  abruptly 
disappears,  resulting  in  a  sudden  increase  or  decrease  in 
propagation  loss. 

Figure  (7.7)  shows  ray  theory  versus  normal  mode  theory 
at  a  depth  of  500  m.  Table  (7.2)  summarizes  the  rays  at  each 
point,  a  similar  pattern  to  that  at  250  m.  Here,  however,  we 
are  out  of  the  caustic  region  at  point  (U)  where  bottom  cut 
off  of  one  arrival  occurs.  And  the  caustic  at  point  (V)  is 
more  obvious.  Figure  (7.8)  shows  the  ray  theory  combined  with 
modified  ray  results,  as  well  as  the  normal  mode  lesults.  We 
see  again  the  interference  of  a  single  ray  arrival  with  a  shadow 
zone  arrival,  point  (W)  to  point  (V).  But  the  normal  mode 
result  shows  an  additional  oscillation,  point  (U)  to  point  (W). 
This  suggests  an  appreciable  shadow  zone  contribution  far  to 
the  left  of  the  caustic  at  point  (V),  much  farther  than 
modified  ray  theory  results  indicate.  This  shows  that  the 
modified  ray  theory  intensity  probably  falls  off  too  quickly, 
and  breaks  down  beyond  some  point  off  the  caustic.  This 
verifies  the  failure  of  modified  ray  theory  far  into  shadow 
zones  discussed  by  Sachs  (49),  The  extra  points  indicated  by 
asterisks  in  Figure  (7.8)  include  in  the  coherent  sum  of  all 
arrivals  the  ray  double  arrivals  to  the  right  of  the  caustic 
at  point  (V),  But  these  tv/o  arrivals  are  in  the  caustic  region, 
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and  the  modified  ray  results  should  be  used  instead.  So  as 
we  indicated  before,  the  modified  ray  results  are  inserted 
in  place  of  the  two  ray  arrivals  in  the  dashed  curve  in 
Figure  (7.8).  This  resuJt  is  clearly  more  appropriate  than 
the  result  indicated  by  the  asterisks  that  includes  the  rays 
from  that  caustic  within  the  caustic  region. 

Figure  (7,9)  shows  the  ray  results  and  normal  mode 
results  for  a  depth  of  1500  m.  Table  (7.3)  summari  zes  the 
ray  arrivals  at  the  various  range  points.  Figure  (7.9)  contains 
a  somewhat  simpler  ray  pattern.  The  caustic  region  and  partial 
double  arrival  region  are  followed  by  a  single  arrival  region, 
another  caustic,  and  finally  at  point  (U),  a  different  double 
arrival  region.  Here  each  arrival  has  passed  through  a  caustic, 
and  one  is  further  surface  reflected.  So  we  have  a  modified 
interference  pattern.  Figure  (7.10)  shows  modified  ray  plus 
ray  theory  and  normal  mode  theory.  Again  modified  ray  theory 
takes  care  of  many  of  the  problems  of  ray  theory.  Beyond 
point  (V)  in  the  normal  mode  solution,  we  see  an  oscillation 
superimposed  on  the  interference  pattern  shown  by  ray  theory. 
This  oscillation  is  similar  to  the  oscillation  in  the  first 
shadow  zone  from  40  Km  to  48  km  and  caused  by  the  matched 
impedance,  discontinuous  gradient  bottom.  We  believe  that  the 
oscillation  beyond  point  (V)  is  also  due  to  bottom  effects. 

Again  the  abrupt  changes  in  propagation  loss  where  bottom  cut 
off  occurs  (points  W  and  X)  are  obviously  regions  where 
improvement  is  needed. 
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We  then  did  the  same  comparisons  Tor  a  50  Hz  source 
frequency.  While  the  ray  arrivals  at  each  point  are  the 
same  (Tables  7.1  -  7.3).  the  frequency  is  different.  So  the 
oscillation  pattern  is  different.  Figures  (7.11  -  7.13) 
show  ray  theory  versus  normal  mode  theory.  Figures  (7.1^  -  7.16) 
are  7e  ray  plus  modified  ray  calculations  for  250  m,  5°0  m, 
and  1500  m.  Also  shown  are  the  normal  mode  results.  For 
250  m,  the  general  agreement  is  good.  However,  the  ray  plus 
modified  ray  result  (Figure  7.1^A)  shows  two  nulls  at  points  (U) 
and  (V),  while  the  normal  mode  result  shows  one  null  near 
Point  (V).  The  nulls  in  the  ray  result  are  due  to  single 
diverging  rays  in  the  caustic  regions  to  the  right  of  caustics 
near  these  points.  We  were  reluctant  to  substitute  the 
modified  ray  result  here,  because  we  didn't  have  two  arrivals  +.0 
remove  as  in  a  normal  caustic  region.  It  is  not  obvious  that 
a  caustic  region  where  bottom  cut-off  of  one  arrival  is 
occurring  is  the  same  as  a  normal  one  with  two  diverging  arrivals. 
It  would  seem  that  this  region  should  be  somewhat  different, 
even  though  the  two  arrivals  couldn’t  be  resolved  in  the 
caustic  region  if  both  were  there.  But  we  put  in  the  modified 
ray  result  in  this  region  to  the  right  of  each  caustic  -  and 
removed  the  single  arrival  from  each  caustic  -  to  see  what  it 
would  do  (Figure  7.1^B).  The  use  of  modified  ray  theory  does 
appear  to  eliminate  most  of  the  first  null,  or  merge  the  two 
together.  It  does  not  work  completely,  which  is  not  unexpected. 

So  our  feeling  that  the  modified  ray  result  is  somev.’hat  like  - 
and  somewhat  different  from  -  a  single  arrival  caustic  region 
seems  to  be  correct. 
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At  500  m  (Figure  ?.15)»  we  again  note  the  general  agreement. 
The  asterisks  starting  at  point  (U)  are  another  set  of  points 
where  we  have  added  the  modified  ray  result  instead  of  the 
single  arrival  in  the  caustic  region.  Here,  however,  the 
region  of  single  arrivals  from  the  caustic  is  quite  wide. 

And  the  modified  ray  is  actually  worse  than  the  use  of  the 
single  ray  associated  with  the  caustic.  This  just  points  out 
that  transitions  from  caustic  to  caustic  at  low  frequencies 
are  quite  complex  and  take  place  over  a  considerable  distance. 

So  caution  must  be  u' ed  in  trying  to  apply  modified  ray 
results  in  wide,  single  arrival  caustic  regions. 

Figure  (7.16)  shows  comparisons  of  the  theories  for  1500  m. 
Beyond  point  (U),  results  from  the  two  theories  start  to  diverge. 
Figure  (7.1?)  shows  calculations  for  the  convergence  zone  and 
adjacent  shadow  zones,  using  normal  mode  theory  and  two  different 
bottoms.  The  use  of  a  hard  bottom  (Figure  7.17B)  not  only 
changes  the  first  shadow  zone  (up  to  point  V)  and  the 
convergence  zone,  but  also  changes  the  region  beyond  point  (U) 
where  the  direct  ray  arrivals  are  getting  weak.  So  we  attribute 
the  disagreement  beyond  point  (U)  in  Figure  (?,16)  to  bottom 
effects  in  a  region  where  non-bottom  reflected  arrivals  are 
contributing  comparatively  less  energy. 

Using  Profile  II,  we  next  tabulated  the  propagation  loss 
partial  sum  as  a  function  of  the  number  of  modes  in  the  sum. 

For  50  Hz  there  are  a  total  of  53  discrete  modes.  For  a  $00  m 
receive!  the  intensity  is  highest,  or  conversely  the  propagation 
loss  is  lowest,  in  the  convergence  zone  ( 50-68  km).  Intuitively, 
one  would  expect  most  of  the  modes  to  contribute  here.  In 
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the  shadow  zone,  the  intensity  is  lowest,  so  intuitively  one 
would  expect  a  fewer  number  of  modes  to  contribute.  So 
from  intuition  alone,  one  would  expect  the  shadow  zone  propagation 
loss  level  to  be  achieved  after  perhaps  10  to  30  modes  were 
summed,  and  the  convergence  zone  level  to  be  achieved  only 
after  almost  all  the  modes  are  summed. 

In  actual  fact,  the  reverse  is  true.  Figures  (7.18)  and 
(7.19)  show  the  propagation  loss  as  a  function  of  the  number 
of  modes  included  in  the  sum.  This  was  done  for  -wo  points 
in  the  convergence  zone  (57.8  and  67  km)  and  two  points  in  the 
shadow  zone  (44  and  50*4  km).  For  both  points  in  the  convergence 
zone,  the  propagation  loss  curve  levels  off  after  about  32  to 
34  inodes  are  included  in  the  sum.  From  this  point  up  to  the 
point  where  the  last  mode  contribution  is  included  in  the  sum, 
there  is  very  little  change  in  the  propagation  loss.  On  the 
other  hand,  the  points  in  the  shadow  zone  show  considerable 
oscillation  in  the  propagation  loss  as  more  modes  are  included 
in  the  sum.  The  propagation  loss  partial  sum  changes  by  as  much 
as  15  dB  from  mode  to  mode,  and  only  with  the  addition  of  the 
last  mode  do  we  reach  the  correct  level.  Thus  the  high 
intensity  levels  in  the  convergence  zone  are  actually  closer 
to  being  average  levels,  and  the  lower  intensity  levels  in 
the  shadow  zone  are  built  up  (or  more  appropriately  broken 
down)  from  the  convergence  zone  levels  by  the  destructive 
interference  of  higher  modes. 
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Table  7.1*  Profile  II,  Receiver  Depth  -  250  m 


Rays  Passing  Through  Range  Points  of  Interest 


Group 

Range 

(km) 

Angle 

Source  Receiver 

History 

I 

59-59.8 

Down 

Up 

Caustic 

Down 

Up 

II 

59.9-63.4 

Down 

Up 

Caustic 

III 

63.6-64.? 

Down 

Up 

Caustic 

Down 

Down 

Far  S.R.*,  Caustic 

IV 

64.4-67.8 

Up 

Up 

Near  S.R. ,  Caustic 

Down 

Up 

Caustic 

Down 

Down 

Far  S.R. ,  Caustic 

V 

68-70 

Up 

Up 

Near  S.R.,  Caustic 

Up 

Down 

Near  S.R. .Caustic .Far  S.R 

Down 

Ud 

Caustic 

Down 

Down 

Caustic,  Far  S.R. 

*  Par  S.R.  -  Surface  reflection  at  convenience  zone  (60-70  Km) 


Hear  S.R.  -  Surface  reflection  near  source  (0-15  Km) 


Caustic  Parameters  for  100  Hz 


Branch/Range 

(km) 

e  e 

0  ( De'g )  r 

AMPDB 

(dB) 

Time 
(Sec ) 

1> 

4- 

A/58.93 

5.919 

82.8 

39.34 

.0025 

0 

B/63.44 

7.427  7.3? 

84 ’.2 

42.32 

.0028 

-7T 

C/64.? 7 

-7.65  -7.55 

84  .  '■ 

»\> 

CO 

.002? 

-7T 

E/79.6 

-I.365  -.66 

88.7 

53 

.00024 

- 
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Figure  7.6  Profile  II,  Normal  Mode  and  Modified  Ray  Theory  (z=250m,  f=100Kz) 
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Figure  7.7  Profile  II,  Normal  Mode  and  Ray  Theory 

0&  =500m,  f=100Hz) 
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Table  7,2»  Profile  II,  Receiver  Depth  -  500  m 


Rava  Passing  Through  Ranee  Points  of  interest 


Group 

Range 

(Km) 

Angle 

Source  Receiver 

History 

I 

56.4-58 

Down 

Up 

Caustic 

•Down 

Up 

II 

58,2-62.2 

Down 

Up 

Caustic 

III 

62.6-65.2 

Up 

Up 

Near  S.R,*,  Caustic 

Down 

Up 

Caustic 

IV 

65.4-69.4 

Ul' 

Up 

Near  S.R.,  Caustic 

Up 

lip 

Caustic 

Down 

Down 

Caustic,  Far  S.R. 

V 

69.6-70 

Up 

Up 

Near  S.R. ,  Caustic 

Up 

Down 

Near  S.R. , Caustic ,Far  S.R, 

Up 

Up 

Caustic 

Down 

Down 

Caustic,  Par  S.R. 

*  Par  S.R.  -  Surface  reflection  at  converKence  zone  (60-70  Km) 


Near  S.R.  -  Surface  reflection  near  source  (0-15  Km) 


Caustic  Parameters  for  100  Hz 


Branch/Range 
( km) 

0  0 
_ °  (Deg)  r  . 

AMPDB 

(dt) 

Time 
( Sec ) 

H 

A/56.33 

4.66 

82.6 

37.63 

.0023 

0 

B/65.26 

7.89 

8.7 

85' 

43.53 

.0029 

-  TT 

C/62,39 

-7.17 

-8.0 

84.7 

41.63 

.0027 

-  TT 

F/8I.05 

.75 

3.76 

92.7 

53.87 

.00027 

-  IT 

_ a _ 
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Table  7.R|  Profile  II,  Receiver  Depth  -  1500  m 


Rays  Passing  Through  Range  Points  of  Interest 


Group 

Range 
( km) 

Angle 

Source  Receiver 

History 

I 

51-53.6 

Down 

Up 

Caustic 

Down 

Up 

II 

53.8-57.6 

Down 

Up 

Caustic 

III 

57.8-58 

Up 

Up 

Near  S.R.* 

Up 

Up 

Near  S.R.,  Caustic 

Up 

Up 

Caustic 

IV 

58.2-69 

Up 

Up 

Near  S.R.,  Caustic 

Up 

Up 

Caustic 

*  Par  S.R.  -  Surface  reflection  at  convergence  zone  (60-70  Km) 


Near  S.R.  -  Surface  reflection  near  source  (0-15  Km) 


Caustic  Parameters  (for-  1.00  Hz) 


Branch/Range 
( km) 

G  6 

0  (Deg)  r 

AMPDB 

(HR) 

Time 

(Sec) 

1 

A/50.95 

3.99 

86 

33.98 

.0019 

0 

C/57.69 

-6.69  -lj.; 

87 

38.42 

.0025 

-IT 

P/87.3 3 

.76  ii,5 

98 

57.98 

.00028 

-IT 

1. 

Figure  7.10  Profile  II,  Normal  Mode  and  Modified  Ray  Theory 

(z  =1500m,  f=100Hz) 
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Figure  7.18  Partial  Sum  of  Modes  (50.4km  and  6?km) 
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Figure  7.19  Partial  Sum  of  Modes  (44km  and  57.8km) 
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8.  Comparisons  for  Horizontal  Caustics 

Up  to  this  point  we  have  concerned  ourselves  with  one 
class  of  caustics.  These  we  have  characterized  as  vertical 
(depth  to  range  slope  IjIOj'S « . 002) .  Now  we  wilj  examine 
the  caustics  we  characterize  as  horizontal  (depth  to  range 
slope  of  1«100  t  t  <  .001).  As  discussed  in  Section  (4), 
these  are  caustics  for  which  a  normal  expansion  off  the 
caustic,  such  as  Ludwig's,  would  seem  more  appropriate. 
However,  for  some  horizontal  caustics,  a  horizontal  expansion 
off  the  caustic  is  still  possible  and  convenient.  For  these, 
we  would  like  to  know  how  accurate  the  shadow  zone  and  caustic 
predictions  are. 

We  originally  intended  to  use  an  arbitrary  velocity 

profile  of  10  to  20  layers  in  treating  this  case.  However, 

the  complicated  ray  patterns  arising  from  such  profiles  tend 

to  interfere  with  the  examination  of  the  horizontal  caustics 

ther  elver.  So  we  selected  a  bilinear  profile* (Figure  8,1), 

for  which  we  could  separate  out  the  caustic  regions  more 

easily.  Figure  (8.2)  shows  the  ray  diagram  for  chis  profile 

and  a  source  depth  of  1000  m.  The  caustics  of  interest  are 

quite  clear.  They  are  summarized  in  Figure  (8.3)  for  the 

second  convergence  zone.  There  is  a  cusped  caustic  (ABC) 

nested  inside  the  caustic  that  limits  the  convergence  zone 

on  the  left  (GH).  Table  (8.1)  sums  up  the  caustic  parameters 

for  caustics  at  two  depths,  201  m  and  281  m. 

#Input  sound  velocity  profile  data  in  Appendix  V. 


123 


NOLTR  7^-95 


At  a  depth  of  201  m,  the  caustics  are  far  enough  apart  so 
that  we  can  examine  the  shadow  zone  for  each  one  individually. 
Looking  at  the  second  convergence  zone,  Figure  (8.4)  shov/s 
the  comparison  between  normal  mode  theory  and  ray  theory  for 
a  50  Hz  source  frequency,  and  Figure  (8.5)  slows  normal  mode 
and  modified  ray  theory.  The  first,  thing  we  notice  in 
Figure  (8.5)  is  that  the  intensity  falls  off  more  slowly 
with  range  into  the  shadow  zone  than  for  the  vertical  caustics 
previously  considered.  For  example,  for  the  caustic  at  153*5  km 
in  Figure  (8.5),  the  propagation  loss  has  increased  10  dB  at 
a  distance  of  3500  m  into  the  shadow  zone.  For  a  vertical 
caustic  (Figure  6.1),  the  first  10  dB  increase  occurs  in  1500  m, 
less  than  half  xhe  distance.  This  is  in  line  with  our  description 
of  horizontal  caustics  as  weak  caustics,  where  the  energy  is 
spread  out  over  a  broader  region  to  either  side  of  it.  For 
both  caustics  in  Figure  (8,5),  the  modified  ray  calculation 
in  the  shadow  zone  is  good  for  about  2000-4000  m  into  the 
shadow  zone.  For  each  caustic,  the  normal  mode  calculations 
show  a  more  rapid  decrease  ir.  intensity  with  increasing  Ar  in 
the  shadow  zone  than  does  the  modified  ray  prediction.  A 
possible  explanation  shows  the  value  of  a  normal  caustic 
expansion.  Consider  point  (M)  in  Figure  (8.3).  We  have 
obtained  a  value  for  the  propagation  loss  there  by  horizontal 
expansion  from  point  (N)  to  the  right.  But  the  caustic  point 
closest  to  point  (M)  is  point  (L^.  and  this  is  a  weaker  caustic 
point  than  point  (N).  So  a  prediction  for  point  (M)  based  on 
point  (L)  would  yield  a  lower  intensity,  one  more  in  line  with 
the  normal  mode  calculation.  We  see  that  one  disadvantage  of  a 
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horizontal  expansion  for  a  horizontal  caustic  is  that  we  may 
be  fairly  close  to  a  much  weaker  caustic  point,  and  not  be 
taking  this  into  account.  But  with  this  caution  in  mind,  we 
have  demonstrated  that  a  smooth,  simple  caustic,  whether 
horizontal  or  vertical,  can  be  treated  by  a  horizontal  expansion. 

Finally,  we  considered  a  deeper  receiver  (281  m),  where 
the  shadow  zone  of  one  caustic  overlaps  the  single  arrival  ray 
region  of  the  previous  caustic.  Figure  (8.6)  shows  a  comparison 
of  ray  theory  to  normal  mode  theory  for  this  depth,  xi,  is 
typical  of  the  previous  cases,  with  divergent  rays  near  caustics 
and  single  arrival  regions  far  enough  to  the  right  of  each 
caustic  (138  krn  to  144  km,  for  example).  We  then  added  in 
modified  ray  results  in  each  shadow  zone  and  caustic  region 
(Figure  8.7).  From  143  km  to  14?  km,  we  see  an  interference 
pattern  resulting  from  a  combination  of  a  real  single  arrival 
and  the  shadow  zone  contribution  from  the  caustic  at  150.5  km. 

The  pattern  is  similar  to  that  in  Figure  (7. 6), around  64  km, 
but  here  the  agreement  is  poor.  The  oscillation  pattern  is 
not  quite  correct,  and  the  level  is  too  high.  This  is  probably 
due  to  an  excessively  large  shadow  zone  contribution  of  the 
type  discussed  in  the  previous  paragraph.  This  dip  •eement 
demonstrates  how  modified  ray  theory  results  using  a  horizontal 
expansion  can  be  incorrect  when  the  point  of  interest  is  far 
enough  from  the  caustic. 
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Figure  8.2B  Ray  Diagram  for  Profile  III .Second  Convergence  Zone 


Figure  8.3  Profile  III  , Caustics  in  Second  Convergence  Zone 
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Table  8.I1  Caustic  Summary »  Profile  III,  Second  Convergence  Zone 


Depth=201  m 


rc 

(m) 

*0 

IBB 

* 

Extra  Phase  Shift 

135455 

10.4 

89.16 

91 

.00111? 

-7T 

153540 

i 

-10.33 

i 

100.85 

93 

.000815 

-3TT/2 

Depth=281  m 


rc 

(m) 

°o 

m 

t 

Extra  Phase  Shift 

133075 

9.9 

87.635 

91.8 

.001105 

-  rr 

150147 

-9.8 

:  98.674 

93.7 

.000798 

-31T/2 
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Figure  8.4  Profile  III,  Normal  Mode  and  Ray  Theory  (z  =201m,f=50Hz) 
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Figure  8.5  Profile  III,  Normal  Mode  and  Modified  Ray  Theory  (z=201m,  f=50Hz) 
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Figure  8.6'  Profile  III,  Normal  Mode  and  Ray  Theory  (z  =28lm,  f=50Hz) 
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9.  Conclusion 

The  comparisons  we  have  made  can  be  evaluated  in  several 
different  ^aysi  in  terms  of  the  usefulness  of  the  modified 
ray  theory  we  have  applied  at  caustics,  in  terms  of  the 
validation  of  a  particular  normal  mode  program  we  have  used, 
or  in  terms  of  the  general  usage  of  ray  theory  and  normal  mode 
theory.  As  in  every  comparison  of  several  theories,  this  has 
been  an  iterative  process.  Knowledge  gained  in  calculating 
ray  theory  curves  leads  to  information  about  normal  mode 
theory,  which  in  turn  leads  to  information  about  modified 
ray  theory,  and  so  on.  In  this  way,  we  feel  that  knowledge 
can  be  gained  about  each  approach  that  would  not  become 
apparent  in  the  consideration  of  one  approach  by  itself. 

It  should  be  noted  that  what  we  call  modified  ray  theory 
(12)  has  also  been  termed  a  caustic  boundary  layer  solution  (13), 
a  uniform  asymptotic  theory  (14),  and  a  caustic  correction. 

To  further  complicate  things,  there  is  a  completely  different 
derivation  termed  modified  ray  theory  (51.  52)  which  is  useful 
when  rays  turn  near  boundaries.  While  this  author  is  comfortable 
with  modified  ray  theory  (it  is  after  all  a  result  of  modifying 
the  basic  ray  equations),  perhaps  "caustic  ray  correction" 
would  distinguish  derivations  intended  for  use  near  caustics 
from  more  general  expressions. 

Whatever  it  is  called,  we  feel  that  these  comparisons 
have  once  more  demonstrated  the  usefulness  of  modified  ray 
theory  near  caustics.  In  many  cases,  the  addition  of  modified 
ray  theory  results  at  caustics  to  simple  ray  theory  yields  a 
much  more  satisfying  propagation  loss  vs.  range  curve  when 
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compared  to  normal  mode  theory.  The  interference  between  "real"' 
rays  and  shadow  zone  contributions  from  adjacent  caustics 
helps  explain  many  of  the  oscillations  present  in  normal  mode 
results  and  absent  in  simple  ray  theory  results.  Thus  caustic 
shadow  zone  contributions  are  an  important  part  of  the 
convergence  zone  intensity  picture.  For  what  we  characterize 
as  vertical  caustics,  the  horizontal  expansion  we  use  works 
quite  well.  For  horizontal  caustics,  the  expansion  can  at 
best  be  said  to  give  only  fair  agreement. 

For  the  horizontal  caustics,  a  normal  expansion  such  as 
Ludwig's  (14)  or  Kratsov's  (15)  should  be  better.  But  in 
reality,  it  should  take  just  a  coordinate  transformation  to 
make  Sachs  and  Silbiger’s  horizontal  expansion  equivalent  to 
Ludwig's  non-uniform  normal  expansion.  The  real  value  of 
Ludwig's  approach  is  in  its  uniform  asymptotic  theory  applications. 
By  uniform  asymptotic  theory,  we  mean  that  for  a  smooth  caustic, 
one  can  obtain  an  expression  that  is  valid  everywhere  -  on 
the  caustic,  in  the  caustic  region,  and  out  into  the  double 
arrival  region.  The  expression  will  automatically  be  identical 
with  ray  theory  in  the  double  arrival  region.  Thus,  haphazard 
agreement  in  the  double  arrival  region  can  be  avoided.  While 
this  approach  is  quite  powerful,  in  practice  it  is  often 
difficult  to  take  advantage  of.  In  many  realistic  deep  water 
profiles  (as  in  our  Profile  II),  the  bottom  cuts  off  one  of 
the  two  caustic-related  arrivals  relatively  near  the  caustic, 
and  the  uniform  asymptotic  theory  is  equally  helpless.  In  our 
comparisons  at  100  Hz,  we  barely  get  one  full  oscillation  in 
the  doubLe  arrival  region  before  bottom  cutoff  of  one  arrival 
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occurs.  For  frequencies  of  1  kHz  or  above  (where  there  are 
many  oscillations  over  the  same  range  increment),  the  uniform 
asymptotic  theory  would  be  useful  as  long  as  the  caustic 
curvature  remained  constant  (otherwise  further  modifications 
would  be  necessary).  However,  for  frequencies  on  the  order 
of  100  Hz,  and  realistic  velocity  profiles,  a  non-uniform 
approach  is  just  as  useful.  This  may  be  either  the  horizontal 
expansion  (12,  13)  or  Ludwig's  (14)  non-uniform  result,  since 
each  has  its  own  advantages.  In  any  event,  the  modified  ray 
theory  we  have  used  does  work  on  the  caustic.  It  works  well 
in  the  caustic  region.  And  for  vertical  caustics,  it  accurately 
predicts  the  near  shadow  zone  field  out  to  a  range  where  the 
pressure  has  dropped  some  40  dB  from  maximum. 

As  far  as  normal  mode  theory  is  concerned,  we  feel 
that  Section  (5)  demonstrates  the  need  for  validation  of  any 
theory  by  comparisons.  Only  this  way  can  apparently  accurate 
calculations  be  verified.  By  comparing  the  normal  mode  results 
to  ray  theory,  we  verified  that  the  caustics  were  where  they 
should  be,  and  that  the  close  in  direct-surface  reflection 
interference  pattern  v/as  accurately  predicted.  So  we  validated 
the  particular  normal  mode  program  being  used.  Comparisons  for 
idealized  profiles  with  exact  solutions  (3,  16)  are  necessary. 
These  standard  cases  give  a  good  indication  of  the  inherent 
accuracy  of  the  program  or  theory  under  consideration.  But 
comparisons  for  realistic  profiles  such  as  we  have  dealt 
with  are  also  necessary  to  verify  a  program  or  theory's 
behavior  for  ranges,  frequencies,  and  cases  of  practical  interest. 


137 


NOLTR  7k-95 


Thus  we  feel  that  an  arbitrary  profile  input  -  such  as  allowed 
in  the  finite  difference  normal  mode  program  we  use  -  is 
desirable  for  evaluation  of  realistic  cases  of  interest. 

One  disadvantage  of  realistic  comparisons  is  that  they 
often  contain  a  bottom,  and  we  are  sometimes  interested  in 
minimizing  bottom  effects  to  examine  other  phenomena,  such  as 
caustic  effects.  For  this  reason,  we  treated  profiles  v/ith  a 
matched  impedance  bottom.  This  way  we  eliminated  first-order 
bottom  reflections.  This  did  lower  the  level  of  bottom  reflected 
energy  in  the  convergence  zone,  and  enabled  us  to  separate 
out  caustic  shadow  zone  effects?  but  we  were  st' 11  getting 
abrupt  changes  in  propagation  loss  where  botton  cutoff  occurred 
and  also  probably  getting  bottom  reflected  energy  from  the 
gradient  mismatch  at  the  bottom.  No  doubt  these  effects  could  be 
accounted  for  by  a  bottom  reflection  treatment  plus  a  diffraction 
correction  past  the  critical  ray,  but  we  were  primarily  interested 
in  caustic  treatments.  In  obtaining  reality  in  our  comparisons, 
we  sacrificed  simplicity.  So  we  were  forced  to  weed  out  the 
phenomena  of  interest  from  other  phenomena  equally  interesting, 
but  not  pertinent  to  this  study. 

With  all  the  problems  associated  with  the  ray  theory 
results?  caustics,  bottom  reflections,  and  diffraction 
corrections  for  bottom  cutoff,  the  question  often  arises  - 
why  bother  at  all?  Normal  mode  theory,  and  other  approaches, 
sum  all  these  effects  purely  and  simply.  A  cynic  might  attribute 
the  continuing  usage  of  ray  theory  to  the  huge  investment  in 
time  and  energy  put  into  ray  tracing  programs  all  over  the 
country.  But  there  are  cases  and  situations  where  ray  theory 
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is  still  better  than  normal  mod"  ;heory.  For  shallow  water 
and  surface  channel  calculations,  normal  mode  theory  is 
generally  accepted  as  superior.  For  deep  waier  profiles,  and 
frequencies  on  the  order  of  100  Hz,  normal  mode  and  ray  theory 
are  roughly  equal  -  even  though  normal  mode  theory  does  include 
all  the  effects  that  must  be  tacked  or  in  ray  theory.  At  higher 
frequencies  ray  effects  become  more  important  than  wave  effects, 
and  ray  theory  is  probably  better.  Ray  theory  is  far  superior 
as  an  interpretive  tool,  indicating  in  ray  diagrams  how  the 
energy  gets  from  one  place  to  another.  And  for  some  applications, 
fast  ray  tracing  programs  (7)  can  easily  outperform  normal 
mode  programs.  Ray  theory  and  modified  ray  theory  are 
especially  useful  in  the  -treatment  of  pulses.  Because  they 
both  have  an  explicit  frequency  dependence,  one  can  trace  ray 
paths  and  caustic  locations  independent  of  frequency.  Then 
one  can  put  in  the  frequency  dependence  as  the  arrivals  are 
coherently  summed  and  the  total  pulse  reconstructed  by  Fourier 
synthesis  (  8  ).  This  is  as  opposed  to  normal  mode  theory  where 
one  has  to  calculate  intensity  for  each  frequency  independently. 
Finally,  the  need  will  always  remain  for  well  known  ray  theory 
calculations  with  which  to  compare  the  new,  more  complex  methods 
of  calculation  continually  being  developed. 
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Appendix  It  Program  RAYTRV  Listing 

(Coherent  Ray  Sorting  Program) 

PROGRAM  RAYTRV 


CDC  6400  FTN  V3.0- 


PROGRAM  RAYTRV { I WPUT»0UTPUT«TAPF2) 

COMPLFX  SUM (1  no  «nu* 

DIMENSION  I RO ( 1 00 i 2) 

DIMENSION  TAPGFT<  110.3) .SP<210»?> .RP(210.2) *VP0(210. 
5  DIMFNSIOM  TEST  ( -T  >  .TEST1  («)  . PUFFER ( 9 . «00 )  .IRFR(9.800) 

EOUI  VALENCE  (8UFFFP.  I.RFR) 


10 


15 


20 


25 


30 


35 


40 


45 


50 


55 


C  TAPF2  IS  OUTPUT  OF  CONGRATS 

C  NCASE  EQUALS  THE  NUMfiER  OF  PROCESS  CAROS  IN  CONGRATS 

C  NCASE  IS  THE  NUMRFR  OF  CASES  TO  BE  PPOCESSEO 

C  IPPN.GT.  0  means  PRINT  OUT  OATA 

C  IPROC.GT.  0  MEANS  AOO  ARPIVALS  COHERENTLY 

C  ISHDAR=N.GT.  0  MEANS  N  SHADOW  ZONE  ARRIVALS 

C  ARE  TO  BE  READ  IN 

CARD  read  ‘STATCMZUTS  -  LtfOK  it,  53 ,24,  9*),  11$" 

READ  976.NCASE 
976  FORMAT (15) 

00  1*20  NCA=1, NCASE 
NSTOP=0 
NSTART=-99 

REAP  Q75. IPRN, IPROC. ISHOAR 
NARPI V=0 

975  FORMAT (315) 

READ  1000. (TEST (J) .J=1.8) 

READ  1000. (TEST1 (J) .J=l.fl) 

1000  FOPMAT(«A10) 

read  (2) (Target ( I ) . i=i »  330) 

READ  (2)  (SP(I) ,1  =  1,420) 

READ  (?) (RP ( I ) ,1=1, 4?0) 

READ  (?)  (VPO(I) .1=1,12*0) 

READ  (?) 

IF (EOF  (?) • ME, 0, )  GO  TO  100 
90  STOP  5 

100  NSTOP=NSTOP+ 1  on 
NSTART=nSTaPT+100 

READ  (?)  (  (BUFFER  (  I  •  J)  .  I=?.9)  .UrNSTART.NSTOP) 

DO  inoi  JKsNSTART.NSTOP 
IF(1RFP(?.JK) .EO.O)  GO  TO  101 

1001  CONTINUE 
IF(MSTOP.NE.ROO)  R 0  TO  100 
STOP  7 

101  NRT=IFIX(TARGET(  110,1)) 

nzt=ifix (Target  (  110.2)) 

NARR I V=  JK-1 

IF (NARR I  V • LE , 0 )  GO  TO  1R20 
I F ( I °PN , L t , 0 )  GO  TO. 725 
PRINT  ?00 

200  FORMAT (1H1 .4nx. lGHCONGPATS  RFSULTS ,////. ?0 X . 

1  34HIMPIJT  INFORMATION  USED  IN  FINDING  ,//.?0  X, 

?  37HRAYS  PASSING  THROUGH  SPECIFIED  TARGET. //.?0X. 

3  17HRANGES  AMO  UfPTHS; 

PRINT  ?  7  5 • T  A  R  G  E  T (  1  OR , 1 ), TARGET (  1  OR «  ? ) . T  A°GFT  (  109.3) 
275  FORMAT (///,SX,?1HS0URCE  RANGE  IN  K YDS= . E 1 5 . A , SX  , 
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PROGRAM  RAYTRV 


CDC  6400  FTN  V3.0-P316  Of 


60 


65 


70 


75 


80 


85 


90 


95 


100 


,  OS 


60 


80 


1  21 HSOURCE  DEPTH  IN  KYDS=.E15.8,5X, 

2  28HS0UHCE  VFLOOITY  IN  KYDS/S=.E15.8) 

PRINT  300,TAPGFT(U7T+?,2)  «T ARGET  (NZT+2.3) 

300  F0RMAT(///,46X.??HSUOf  6CF  DEPTH  IN  <Y()S=  *  E 15  .  ft  *  5X , 

1  27HSURF ACE  VELOCITY  IN  K YOS/S= , E?0 . 8 ) 

PRINT  380 » TARGET (M7T+ 1.2) «  TARGET (N7T+1 ,3) 

350  F0»FAT(///,46X,21HPOTT0m  DEPTH  IN  KYl).S=,E15.8,5X, 

1  26H80TT0M  VELOCITY  IN  K YDS/S= . E20 . 8 ) 

IFITAWGETI  108.3) .LF.o.)  GO  TO  3690 
PRINT  375. TARGET!  1 08 . ?) , T ARGFT  (  108.1) 

375  FORMAT <///.36X. 178FPF0t)ENCY (RAO/S) =.E20,8.5X, 

1  22HATTENUATION  IN  OR/K YO= . E20 . « ) 

3690  NVP0=IFIX (VPO (210. 1 ) ) 

3700  PRINT  370? 

3702  FO°MAT (///.45X.28HCOMGPATS  VFLOCITY  PARAMETERS. 

1  /  A5X  28H - - 

2  //.6H  LAYER. 7X.AHZ0-KYD, 12X .8HV0-KYD/S. 14X.2HG0, 17X.2HG1 ♦ 

4 1 7  X .  2HG2, 14X.7HV-KYD/S//) 

PRINT  3712.  (N.  (VPO  (M.  J)  ,  J=1 ,6)  .N-1'.NVPO) 

3712  FORMAT  (14. 1R6"F1  9.8) 

PRINT  400 

400  FORMAT*///, SO*. 21HTARGET  DEPTHS  IN  KYDS) 

PRIM)  8')0  .  (TARGET  (I,  ?)  ,I  =  1.NZT) 

500  FORMAT  .  i  )*  ,  6E20 , 8 ) 

PRINT  550 

550  FORMAT (///.SOX, R1HTAP6RT  RANGES  IN  KYDS) 

PRINT  500.  (TARGET (1,1) ,1  =  1, NR T) 

PRINT  700 

700  FORMAT*/////,  4nx,?2HC0LLECTI0N  OF  ARRIVALS) 

725  DO  1350  J=1,NZT 

1350  TAPGFT  ( J,2)  =TA°r,FT  (.1,?) /TARGET  *  104. cj 
00  1360  J=l,NRT 

1360  TAPG^T(J,1)=TAR0PT(J,1)/TARGET(  104,1) 

DO  1370  J= 1 , N APR T V 

RUFFFR(3.J)=RijFFFP*3,  J)  /  ( 1 . 74532925E-2 ) 

1370  RllFFFP  (4.  J)  =  Hi  jFFrR(4*J)/(  1,745329255-2) 

IF  ,  I SHOAR.LE . 0)  GO  TO  740 


C  ALL  QUANT ITIES  READ  IN  HERE  ARF  iN  ORIGINAL  UNITS 
C  SAME  UNITS  AS  FFED  INTO  CONGRATS 

C  ZC  IS  DFPTn  OF  INTEREST 

C  RC  IS  CAUSTIC  manor  AT  DEPTH  ZC 

100  £  (HUHnetiiactumotnnjMtiif.siioinHuHms 

C  MAKE  Zf.P C  HOT  H  SAMP  UNITS 
C  EITHER  METERS, FEET, YDS, FTC 

Q  aoaaaaaaaaafiaacaaaaaaaoaaoaaaaaoaa 

C  THFATO  IS  INITIAL  A  NGl.P  ( I M  OFG)  OF  RAY  GOING  THPOUGH  CAUSTIC 

C  TC  IS  TRAVEL  TIME  TO  CAUSTIC 

C  APPRO  is  PROP.  LOSS  TO  CAUSTIC  IN  OP 

C  PS  IS  ANY  ADDITIONAL  PHASE  SHIFT-DUE  TO  SURF/.CF.ETC, 

C  W1  IS  THIRD  D PRIVATIVE  OF  w  NF.PDED  IN  AIRY  FUNCTION 

C  ISIGN  IS  SIGN  OF  THIRD  Of R I V^T I VE--* 1  OR  -1 

l1*  5 


i 


‘5 

| 


no 
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PROGRAM  RAYTRV  CDC  6400  FTN  V3.0-P316  OP 

C2=TApGFT(10O,3)/TARGFT(104,l) 

AK=TA«GFT  ( 106,3>/C2 
AIO=A I (0.) 

DO  2990  IF=1,ISH0AR 

HR  H5  RF  AD  201  0 , THF A TO. 70  «RC*TC*AMPP8,wi  » PS »  ISI6N 

2010  FORMAT  (7F10.3,I?) 

XIC=COS(T>'FATnM  .7AS32R25H-2) 

00  2R00  IF  2=1 • N7T 

IF'ARS(TARGFT<IF2.?)_ZC). PE. .00001)  GO  TO  2900 
120  00  2«P9  IF1=1.nRT 

DELR=T  ARPET ( IF l *  1 ) -RO 

RHO=  (  { .66*7)  )  )  tffiFLR0 1  SIGN 

IF  (RHO.l.E  .  0,  )  GO  TO  2R99 

N  A  R  P  T  V  •"  (“  a  R  R I V  ♦  1 

125  IF(NAR:mv.LF.AO0)  GO  TO  2200 

STOP  4 

2200  RDFFF" (PiAJARPJV) =AmpdR-20.*ALOG1O (A8S (AI (RHO) )/AlO) 

8UFFFP  (S.MAPRIV)  =TOXTC*QELP/C2 

BUFFPWI?  ,NA;.RI  V)  =TA»r,FT  ( 10  8.3)  *8UFFFR  (5*NARRIV)  -. 785398 16+P‘ 
130  IRFR (ft.NARpl V) =1^1 

I R  c  P  (2«NARRIV)  =  I  F  2 
BOFFFP ( 3 « NARRI V) =n . 

BUFFPR (A.NARPIv) =0. 

IRFP (O.MAPRI V) =9PPQP 
13S  ?R99  CONTI NUF 

2900  CONTINUE 

2990  CONTIMUF 


140  140 

740  00  750  I=1.MARRIV 

750  IRFR(1,I)=100»IRPP(2»I)+IPFP(8,I) 

J=9*MARRJV 

CALL  COMSOT(IRFR.J,o.O,l,l) 

145  IF (IPPM.LE.O)  GO  TO  1400 

POINT  1200 

1200  FOo^AT(1H-,3X.1?HTARRET  DEPTH, 15H  TARGET  PANGF  , 

1  15H  INITIAL  ANDLF  , 1 5H  FINAL  ANGLE  ,15H  TRAVEL:  TIME  » 

2  1RH  PROP.  LOSS  *  1 5H  PHASE 

1 RO  3  1 RH  PHASE  SHIFT  tlOHNO.  OF  PFV) 

POINT  1250, (TFST(J) ,3=1,4) , (TFST1 (L> *L=1 .4) 

12*0  FORMAT (1H  « 1 2 A  1 0 ) 

DO  1390  J=1,NARRJV 

PS=RUFFF.P(7,  J)-TARRf  T(  10R,3)  »RUFFFR  (5,  J) 

1 S5  1390  PPINT  139G,TARG?T (IRFH (2, J) ,?), TARGET (IHFR(3,J) *1) , 

1  .(PUFFFR  (I,J).I  =  3.7).PS,HFR(0,J) 

1395  Format  (RF.  15.  R.  II  0) 

1400  I F  ( I P°OC . LF . 0 )  GO  TO  1R20 
DO  14?0  0=1.100 

160  1420  S'JM(J)=CMPLX(0.,0.) 

J-0 
1C  =  0 
IZT  =  0 

1 4R5  IZT  =  IZTH 
165  lo5  1460  I RT= 1 
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PROGRAM  RAYTRV  COC  MOO  FTN  V3.0-P316  OPT  = 

1465  IC=I01 
1475  J=J*1 

IF (J.GT.MARPIV)  GO  TO  1700 
lF(I°FR(2»J).Mr.I7T)  GO  TO  1*00 
170  ±?0  IF(IPFR(ft»J) .NF.fRT)  GO  10  1590 

C - THIS  FLAGS  POTTO**  (  OP  SURFACF)  POUNCES  VMEN  LOSS 

C - - - ! M  CONGRATS  IS  SFT  TO  9.E+03  OP  SO  ONE 

175  C - CAN  DELETE  THEM  FROM  THE  SUM  OF  ARRIVALS 


180 


IPS 


1477  CONTINUE 

IF (PUFFER (6*  J) .GE.9.E+03)  GO  TO  1485 


190 


195 


200 


205 


210 


21P 


220 


195 


220 


AMP=10.*» (- (PUFFFR (a, j)  /?0. )  ) 

PH=  PUFFER (7 • J) 

Rr=AMR*COS (PH) 

AIG=AMP«SlN(PH) 

DIJM=CMPLX  (PE.  AIG) 

GO  TO  14R0 

14«5  OUM=CMPLX (  0  . .  0  .  ) 

1490  SUM(IC)=SUM(IC)*OUM 
IRO ( IC* 1 ) =17 i 
IRD ( IC»  2) =IRT 
GO  T0  1475 

1590  IF(IPT.EO.NRT)  GO  TO  1455 
IftT=IRT  + 1 
J=J-1 

GO  TO  U65 
1600  J=J-1 

GO  TO  1455 
1700  PR  TNT  1750 

PRINT  17*>0  .  T^ST  ( 1)  .TEST  (2)  .TEST  (3) 

00  1°00  IJsl.IC 

IF ( (REAL (SUM ( I J) ) .FO.O. ) . AND. ( A  I M AG ( SUM ( I J ) ) .EO.O. ) ) 
2100  FOWr- A T  (110«4E15.7) 

Y  =  CAOS(SU'MIJ)  ) 

Y1=PFAL(SUm(IJ) ) 

Y?  =  A  IMAM  (S'l*  (  I  J)  ) 

P1=ATAM2(Y2.Y1 ) 

OR=-?0.‘*ALOG10  (Y) 

GO  TO  1799 

1750  FOC‘MAT(////*30X.27HPFSULTANT  aT  TARGET  POINTS.//* 

1  3*  < 1 ?hT  ARGE T  ofRTH.AX. 15HTARGET  RANGE  .6X» 

p  ]  rurfruLT anT  A'-plITUOf,  iox.hhphaSF.) 

1740  Format  (/.3A  10.7 x.3-*9SI .  12Xi2HOP«  12X«3HRAO) 

1770  Y  =  oqq^vj. 

OP-99999. 

1^7 


GO  TO  1«0 
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PROGRAM 


22*  c25 


RAYTRV 


CJ>C  6400  FTN  V3.0-P31b 


P1=0. 

1799  PRINT  1  RIO,  TARGET  ( IPO  ( I J,  1 )  *2)  ,  TARGET  ( IRO  ( I-Jt  21  ,1)«Y»D8.P 

1800  CONTINUE 

1810  FORMAT (/,5E15. 7) 

1820  CONTINUE 
ENO 


148 


S0ouce  SaNGE  In  KyOS*  0.  SOuRCE  DEPTH  In  KyDS«  . 33333?67E»00  ROijRCE  VELOC  ;Ty  IN  KyDS/S"  ,166500l0E*0I 


•“•'-^r'oooc  o 

c*  .i  c-  o  c  w  r*  r 

|  |  |  4  4  4  ♦  ♦  ♦ 

U  a)  . i*  *.  jJ  u*  a.  a.' 

O  n  *1  -"l  > 

>-  O'O  jnno 

*■  >-^^OiPr'fv:CKa' 

cj  O'^nno'^occ' 

r\  a  r  r  r  '»•  ^  ^  si  «»> 

o  cr  n  n  ^  o  c*1 
?<*r-rv^— 


*i  cvi  fv>  <\i  ru  cv»  cnj  <\j  f\/  r\i 

♦ 

U  UJU'LJUJUJUJUJUJUJ 

h-  K  '•  "IN  erf-.  -  — 

•c 

•*«  «ccr>c  r*) 

<•  m  *  *  —  r  v  a  •  r  i 

o  >nn>f  n  ^ 

«»  r  *>  r.  - 

•o 


« 

t 


1U9 


COLLECTION  OF  ARRIVALS 

Arrivals  with  99999  in  Last  Column  Are  Shadow  Zone  Arrivals  from  Appropriate  Caustics 

TARGET  r.T3TH  .TARGET  RANGE  INITIAL  ANGLE  FINAL  ANGLE  TRAVrL  TIME  PROP,  LOSS  .PHASE  PHASE  SHIFT  *N0.  OF  REV 

u  KM  .  deg  PEG  Sec  08  '  RAO  RAO 

Soo-ico .«>POOOOOOE*nS  -.9i<S332mE»fll  .lASaG^SaEtoa  .«032?T02E*02  .Vo^SOTSoEtOA  .126Gl‘>G7E*-S  -.G28310S3E*O1  3 
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NNW^nNNWonNNfyo'nNWwcr  m  rg  c\j  c\j  o  m  ru  nj  ru  o  nAjcuc'Jorjrurjfgo  mr\jrur\jor)f>j<\jfwo  t 
0.0  o  o  o  o  o  o  o  o 

OOOOOOOOOO 

o  o.o  O  O  O  0.0  o  o 

OOOOOOOOOO 


OOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOOO* 
(JUJlijUJU.1  UU’U.UJUU.  UJUJUJ(l.UjU.UJU,‘UJUJLiIU.  LkJUJLuU-'  LiJUJLiJUJU.  U  A'  U  UJ  U.'  U  UUd  liJ  U1  U-U-Ui  U.  LJ  UJ  U-  . 

or)h*oricm^^.r(cr>h»cfr)on^^,fncrih-^>r>cnr«.'£r^or;r»<r>cr. 

<CO  O  ©  '3  CC  O  O  X  X  X  O  O  X©XOO't-©XOOX'<OT*00.»©XOOCt>X>CGOOXT'<COO  C  X>  COO'S- 

ns(AO,«r>N;AO'"*nNu,i7.«r>Mfl(>^oNino— 

(Me— —  —  o^nnroc^nr,  rjc^nnwo^mnAc-nnPjc^ni 

^  ih  s  »-n  -f  ia  ®  k  ■#  3  »-•  n  ^  i/i  (O  r*N  fi  x)  — <  r—  *7  ^  s  4ioo  *-•  f-  «»  .n  ©  *“•  t** 

rH  o  fu  n  i/i  *-•  aj  ruh*  in  «-«  x>  rvif^in^o  wn  in  — ♦  o>  rv>  in  — •  x.  ru  r-  in  •~*©njr-in«-«®(\jf'-in«-*©i\jf'*in  i 

■f  ^FJN'O  ^"1  N  o  >0  4  -»nNO^-«r)N'0-tHnN^ 


»  »  »  •  »  i  i  »  *  i  i  :  i  i  i  i  I  i  i  i  »  i  i  i  i  »  •  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i  i 

ininininintninininininininmtninininininintnintninintninininininintnintninintninininintntnininintni 

*’  i  i  .  '  c  :  -  (  (  *  c  v  c  c  l  :  :  c  '  :  L  t  i  •  f  -  r  -  i.  (  c  l  :  .  v  T  *  i  <  ^  r  .  c  « 

♦  ♦  ♦  ♦  ♦  ♦♦♦♦♦♦  ♦  ♦♦«♦♦♦*♦  +  ♦»  +  « 
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Appendix  II «  Normal  Mode  Program  Listing 


The  following  is  a  card  by  card  listing  of  the  data  deck. 
This  precedes  the  program  listing.  While  the  program  has 
undergone  extensive  revision,  the  main  logic  and  substance  of 
the  program  follows  Newman  and  Ingenito  (22). 


Card  1  READ  500,  KKK,  INC 
500  FORMAT  (214) 

KKK  -  number  of  sound  speed  profiles  for  which 
calculations  are  to  be  made 
INC-If  M0DSHAPE=1,  -prints  amplitude  at  layers: 

1,  1  +  INC,  1  +'2  •  INC,  etc. 

Card  2  READ  4500,  PROPLT,  MODSHAPE,  MODEPLT,  GRUVEL,  GRUPLT, 
GRUVELS 

4500  FORMAT  (F5.1,  215,  3F5.1) 

PROPLT  =  1  -  plots  sound  speed  profile  otherwise  - 
no  plot 

MODSHAPE  =  1  -  prints  mod  amplitude  values 
otherwise  -  no  print 

MODEPLT  =  1  -  plots  mode  amplitudes  -  otherwise  - 
no  plot 

GRUVEL  =  1  -  indicates  group  velocities  are  to  be 

calculated  otherwise  -  program  terminates 
GRUPLT  =  1  -  plots  group  velocity  curves  otherwise  - 
no  plot 

GRUVELS  =  1  -  prints  group  velocities  over  frequency 
range  specified  otherwise  -  no  print 

Card  3  READ  1000,  TITLE 

1000  FORMAT  (10AS) 

TITLE  -  word  and/or  number  description  of  problem 

Card  4  READ  2000,  VELNO 
2000  FORMAT  (I 10) 

VELNO  -  number  identification  of  sound  speed  profile 

Card  5  READ  4000,  F,  CT,  Rl,  R2,  H,  LI,  ND,  NM,  NF,  EPSILON 
4000  FORMAT  (5FIO.3,  4l4,  F10.4) 

F  -  Frequency  of  source 

CT  -  second  layer  sound  speed 

Rl  -  density  of  first  layer 

R2  -  density  of  second  lay°r 

H  -  fi^st  layer  de^b  (water  layer) 

LI  -  number  of  .increments  into  which  first  layer 

is  to  be  divided  for  finite  difference  equations 
ND  -  number  of  sound  speed  depths 

NM  -  number  of  normal  modes  or  eigenfunctions  wanted 
NF  -  number  of  additional  frequencies 
EPSILON  -  criterion  for  acceptable  solutions 
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Card  6  READ  6000,  (Z1(I),  C1(I),  1=1,  ND) 

6000  FORMAT  (2F10.3) 

Z1(I),  Cl( I )  -  profile  depth  and  sound  speed 

Card  ND+6  READ  7500,  SD,  (RD(I),  1=1,5) 

7500  FORMAT  6F10.3 
SD  -  Source  Depth 
RD(I)  -  Receiver  Depth 

Card  ND+7  READ  14000  (LB1(I),  LB2(I),  1=1,10) 

14000  FORMAT  (10(?.l4)) 

LB1,  LB2  Specific  modes  to  be  calculated 
i.e.  1  5  100  110  0  0 

Means  calculate  inodes  1  through  5,  and  100 
through  110 

Card  ND+8  READ  1B00  F 

1300  FORMAT  (F10.3) 

F  -  New  source  frequency 

Card  ND*NF+8  READ  14000,  FMIN,  FMAX ,  FDELF 

FMIN  Lowest  frequency  for  group  velocity  calculation 
FMAX  Largest  frequency  for  group  velocity  calculation 
DELF  Frequency  increment  for  group  velocity 
calculation 
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PROGRAM  NQRM003  COC  6400  FTN  V3.0-P316  OPT  =  l 

PROGRAM  NORMOOK  InpijT.0UTPUT*T4PF1*TAPF7,PUNCH=TAPE7,TAPF9) 
COmmom/1/  XT (  3100) 

CO,4MON'/?/Z2M  (3100)  ♦  C2D  (3100) 

COMmom/aa/CP (20,20) , V ( 20 , 20 ) ,6V ( 20 , 20 ) »C ( 20 , 20 ) , DKN (20 , 20 ) 

5  COMMON/ A-T/Z 1  (50)  ,C1  (50)  »  Z  IN  (50  )  «LL*(10)  «LUH(10)  ,MMU<30)  ,LR1  (10)  , 

*LR?(10) 

COMMON/AC/AKN (  12) «RR (  12),PV(  12)»AN(  12),7F<12) 

COMMON/ AO/EPS  I  LON, H, HD. HD  3 *OL ,DLS, R , I T , CT, CTD. MOOSHAP  ,M0DEPLT, 
9GRUVEL.GRUPLT.  IZEPO,  jrv 

■)  COMMON/AE/LOEPTH  .LAMP  ,lgru  »LFRE 

DIMENSION  INDEX (6) , AMPMOO ( ft . 12) 

DIMENSION  RD(S) 

DIMENSION  ZT (3100) ,ZZ (3100) 

DIMENSION  PLTmOOE ( 264 ) .TITLE(IO) 
r  EOU7VALEMCE (??n,ZT) , (XT.ZZ) 

data  LAmP/QhamplITDDE/ 

DATA  LOEPTH/SH  DEPTH  / 

DATA  LFRE/9HFRE0UFNCY/ 

DATA  LGRU/RHGROUP  V EL/ 
n  11111  FORmaT(?I4) 

000  F0PMAT(F5,1.?I5.3F5.1) 

1000  FORMAT ( 1 0 A3 ) 

2000  FORMAT ( 1H1 . 10AR) 

3000  FORMAT (110) 

S  4000  FORMAT!///) 

5000  FORMATS  FREOUENCY/BOTTOM  V./H?0  DEN. /BOTTOM  DEN/H20  DEPTH/  LI/  ' 
«/  nm/  ne /  EPSILON  o) 

6000  FORMAT (5F10. 3,414, F10. 4) 

7000  FORMAT (2F10. 3) 

0  30  P000  format (9  SOUND  SPEED  PROFILE  *> 

«000  FORMAT ( *  DEPTm  VFLOCITY  9) 

10000  FORMAT (1H1, 9  SOURCE  FREQUENCY  =  *,FB.3> 

11000  FORMAT)/) 

12000  F Odm at ( *  MAXIMUM  ND.  OF  ZERO  CROSSINGS  =  *,I4) 

5  13000  FO°maT ( *  mode  CUTOFF  AT  THIS  SOURCE  FREQUENCY  «) 

14000  format (10(214)) 

15000  FORMAT!®  MODE  =  »,I3) 

16000  FORMAT ( *  WAVE  NUMRER  PHASE  VELOCITY  *  R  *) 

17000  FORMAT (3E20. 13) 

0  1P000  FORMAT ( 1 H 1 ) 

19000  FORMAT (*  MODE  AMPLITUDES  FOR  SOURCE  FREQUENCY ( HZ )  =  *,FA.2> 

200  0  0  FORMAT (10X, 12 (3X,9modf»,  13) ) 

21000  F0PMAT(4X,A.  ,12(X,A9) 1 
22000  FORMAT ( 1 3F 1 0 . 3 ) 

5  23000  FORMAT (F10. 3) 

24000  FORMAT (3F10. 3) 

25000  F0RmaT(«  GROUP  VELOCITIES  *> 

26000  FORMAT ( X, A8. A1 , 12 (X.Ap.Al ) ) 

27000  FORMAT (ion  1 .4) 

.0  27001  FOPMAT (2F1 1 .4) 

CALL  PLOTS (RLTmodE, 254, 1 ) 

C  FIVE  RF.An  STATEMENTS  STORF  CONSTANTS  AND  FLAG  INFORMATION 
RE  AO  1 1 ) 1 1 , KKK , INC 
00  9997  M=1 , KKK 

>S  55  READ  flOO.PROPLT.MODSHAP  .MODF.PLT  ,  GRUVEL  *  GHUPlT  ,  GRUVELS 
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,0 


iS 


'0 


15 


<0 


^5 


30 


Qt; 


00 


n? 


[0 


READ  1  000 « TI TLE 
PRINT  2000. TITLE 
READ  3000  *  V£LNO 
IF<PROPLT. £0.0.0)  GO  TO  5 
60  CALL  PLOT (0.*9.*3*1) 

CALL  SYMHL4 (0.0, -A. 0,0, 14* 10HPROFILE  10,90.0,10) 

CALL  NUM«R(0.0, -4.0,0. 14, VELNO»®0.0»-1  ) 

CALL  PLOT<2.6,0.*3,1) 

5  PRINT  4000 

PE  AO  6000 *F,CT*R1*P2,H,LI. NO, NM,UF, EPSILON 
PRINT  5000 

PRINT  6000  * F, CT, R1 ,P2*H,L I,ND,NM.NF, EPSILON 
PRINT  4000 

C  SIXTH  READ  STATEMENT  STORES  SOUND  SPEED  PROFILE 
READ  7000* <Z1 (I) *C1 (I) ,I=1*ND) 

PRINT  8000 
PRINT  9000 

PPIMT  7000* (Z1 (I) ,C1 (I) ,I=1,ND) 

READ  7500*SO, (°0 ( IL) *IL=1,5) 

7500  FORMAT (6F10.3) 

ISTEP=H*INC 
JSTEP=5»INC 
SDEP=SO/H 
DO  7  1=1*5 
7  RO(n=RD(I)/H 

PUNCH  7600 ♦ SOEa  * !PD ( IL ) . IL=1 *5) *LI 
WETFST=.6/FL0AT(LI) 

7600  FORMAT (6F10.6* 15) 

C  NORMALIZATION  OF  SOUND  VELOCITY  PROFILE  DEPTHS 
00  10  1  =  1, NO 
10  Z1MI)=Z1  (I)/H 
N=L I + I 
0L= l . 0/L I 

c  normalization  OF  incremental  DEPTHS 
DO  20  1  =  1  * N 
20  Z2N(I)=DL*(I-1) 

C  CALL  TO  SUPROUTIMF  '.HICH  LINEAR  INTERPOLATES  BETWEEN  THE  ABOVE  TWO  SI 
c  NORMALIZED  rater  depths  TO  yield  A  set  OF  SOUND  VELOCITIES  FOR  Ti 
C  NORMALIZED  INCREMENTAL  WATER  DEPTHS 
CALL  ZCINTEK  <ND*N.F,C2MIN,PR0PLT) 

1=0 

OLS=DL«Dl 
CTD=CT«CT 
HD=H(tH 
R=R2/R 1 
IT=N 

HD3=HD6H 

C  CALCULATION  OP  LARGEST  KN  { AK? )  AND  SMALLEST  KN  ( AK1 ) 

30  AK1=(6.2R31B53071 A«F)/CT 

AK2=(6.2B31B63071°*F) /C2MIN 
FO=(3°.47B41760436*F«F) 

AK3= AK 1 

■IF  (MODEPLT.EO.O)  GO  TD  4C 
CALL  PLOT (5,*0,*3*1) 

CALL  SYMriL4(0. 0,-6, 0,0. 1 4 , 20HSOURCE-  FREQUENCY ( HZ ) *90.0*20) 
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5 


?0 


*0 


55 


o 


5 


o 


5 


r* 


5 


CALL  NUMBRtO.O. -P.0.0. 14. F,e>0.0«l.  ) 

C  CALCULATION  OF  MAXIMUM  HOOF  AVAILABLE  USING  A  SUBROUTINE  WHICH  ITERA 
•  C  FROM  POTTOM  TO  SURFACE  USING  A  FALSE  POSITION  TECHNIQUE 

C  SUBROUTINE  ITEPATE  CALCULATES  A  MODF  SHAPE  FOP  A  GIVEN  EIGENVALUE 
115  C  BY  THE  METHOD  OF  FALSE  POSI TION (REGULA  FALSI) 

40  IZF.RO=0 

CALL  ITERATE < AK1 ,F0,H,C2D»DLS*HD»HD3»CTD*lT»NCR*R»Z»A,ZZ»0L»IZERf 
*1) 

DO  41  IL=lt6 

41  I Mf>E X  ( IL )  =  0 
PRINT  10000, F 
PRINT  11000 
PRINT  12000, NCR 
PUNCH  11500, F,R1,H, NCR 

11500  FORMAT  (3F10.3, 14) 

PRINT  11000 
ir(WCR.EG.O)  50,60 
50  PRINT  13000 

READ  14000, (L91 <I) «L8? (I) ,1=1,10) 

GO  ro  250 
60  !COUNT=JPAGE=0 

READ  14000,  (LB1 ( I ) ,  Lp2 (I), 1  =  1, 10) 

PRINT  14000, (LB1 (I) ,L«2( I) ,1=1,10) 

DO  70  1=1,10 
LLP ( I >  =Lftl ( I ) 

70  LU«(I)=LB2(I) 

FLAG=0 . 0 


145 


SKIP=IC=1 
DO  120  J=l,10 
JJ=LLP ( J> 

AO  IF ( JJ.LF .NCR)  110,90 

«0  LUP(J)=NCR 
GO  TO  120 

110  IF ( JJ.EQ.LU0 ( J) )  120,115 

115  JJ=JJ+1 
GO  TO  80 
120  CONTINUE 
130  I  COUNT  =ICOUNT  +  l 
AK1=AK3 
I=LLH( ICOUNT) 

IUM=LUR ( I  COUNT ) 

I F  (  I  .EQ,0''>  240*140 

140  PRINT  15000,1 

C  SUAROUTINF  HALF  DETERMINES  EIGENVALUE  FOR  EACH  MODE  SHAPE  BY  THE 
C  HALF  INTERVAL  SEARCH  TECHNIQUE 

CALL  HALFIAKl » A*? ,FQ t NCR , A, I , IC> 

RA=SORT (A) 


LJ=LI~1 

0OD=EVEN=0.0 

C  CALCULATION  OF  NORMALIZATION  CONSTANT  USING  SIMPSONS  RULE 
DO  160  J=2 , L I , 2 
160  EVEm=EVFw+ZZ(J)*Z7(vIS 
DO  170  J=3»LJ,2 
170  ODD=ODD+ZZ (J) «7Z (J) 

165  C  SIMPSONS  FOLIATION  PARABOLIC  PULE  VERSION  HILDEBRANO 
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:0 
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'0 


15 


0 


lc 
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N0RM0D3  cnc  6400  FTM  V3.0-P31A>  OPT  =  l 

AIl=R]°'(0L/3.)<>(ZZ(l),:,ZZ(l)+4.*EVEN  +  ?.'i,00D  +  ZZ{N)',>ZZ(N)) 
A1?=P?#(1.0/(2.0*RA) ) 

AN0=(1.P/(AI1*4I?> ) 

AN(TC)=SORT (AND) . 

C  SUHROUTINE  PLOTS  MODE  SHAPES 
CALL  SHAPE (N, IC* I ) 

RR  (TC)  =  i6.28318S30'5F'&P2<>AMD)/(CT<>AKN(  IC)  *2.0*RA> 

PV (IC)  =  <6.2831H530°F>/AKN( IC) 

PRINT  18000 

PRINT  1700  0  « AKN  ( IC )  .PV(IC)  ♦RP(IC) 

PUMCh  1 700 0  *  AKN ( I C ) ,PV(IC) *RR(IC) 

PRINT  4000 
JPAGF=JPA6E+1 
IF ( JPAGE.EQ.7)  180*185 
180  PRINT  18000 
JPAGE=0 

C  TEEf  TO  DETERMINE  WETHER  ADDITIONAL  MODES  ARE  REQUIRED 
188  AK?=AKN(IC) 

AK1=AK3  ' 

MNU  ( I C )  —  I 
IC=IC+1 
1  =  1  +  1 
60  TO  230 
190  PRINT  18000 

IF  (  (AFLAG.EO. 1 . ) . and. (FLAG.E0.2. ) >  GO  TO  250 
AFLA6=1. 

ENOFILE  7 
JP4GE=0 
I K= I C— 1 
IC=1 

222  JPUN=0 

PUNCH  21400,ZT  (INDEX (1) } ♦ < AMPMOD ( 1 ♦ K2 ) *K2=1 *12) * JPUN 
228  JPUN- 1 

DO  1 1 RO  111=2*1 NDEN 

11Q0  PUNCH  ?)4 00* ZT( INDEX (III) >  ,  OAMPMOOUI  I.K3)  *K3  =  1 » 12)  » JPUN 
21400  FORMAT (F 1 0 »6»  3F20 . 1 3/4F20 . 1 3/4E20 • 13/E20 •  1 3 » 1 2 ) 

ENnFILF  7 
PRINT  18000 

IF  (FLAG. ED. 2)  GO  TO  250 
GO  TO  236 

230  1F(I.LF.(LLP<1)+1)>]?30,1300 

1230  00  1235  * N 

IF  (  (ZT (J)-WFTFST) .LT .SDEP. AND. (ZT ( J) +WETEST) .GE.SDEP)  INDEX ( 1 ) =J 
DO  1?3P  I LK= 1 .5 

IF ((ZT ( J) -WETEST) .LT.PN (ILK) . AMD.  ( ZT ( J) + WETEST ) .GE .RD ( ILK) . AND. f 
«  (  ILK)  .NF.O. )  I MHEX ( ILK* 1 ) = J 
1232  CONTINUE 
1235  CONTINUE 
I NOFN= 1 

DO  1236  ILKl-2»6 

IF  ( INDEX  ( ILK  1  )  .NE.O)  INI)EN=INDEN+1 
1236  CONTINUE 

1300  AMPMODU.  (IC-1)  )  =  X  T  (I  NOEX  ( 1 )  ) 

DO  1303  JJ=1 .N.ISTEP 
IENO=JJ* JSTEP 
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1303 


225 

130b 

236 

240 

250 


260 

270 

?«0 


2R0 


2^5  300 
310 
320 


330 


34  0 
350 

C 

360 


370 


3«0 


275 


IFN0=MIN0(IEN0,N) 

PRINT  27000*  <ZT(LL» , XT <LL) ,LL=JJ, IEND, INC) 

CONTINUE 

aflag=o. 

DO  1305  1 1 1  =  2, INOEN 

AMPMOiXIII, (IC-1) )=XT( INDEX (III)  ) 

IF (IC.P0.13)  DO  TO  1 RO 
IF(I.F0.IUP*1)  130,140 
FL AG=2 
GO  TO  1R0 
NF=NF-1 
FNDFILE  7 

I F (  NF  .GE.O)  260,270 
RpAO  23000  »F 
GO  TO  30 

IF(GRl)VPL.EO.l.O)  200,9997 
READ  24000  ,FMIN,FMAX,DELF 
IC=1 

K=( ( <FNAX-FmIN)/OELF) *0.5) +5.0 
I  COUNT  =  J=0 
I  COUNT  = I COUNT ♦ 1 
I=LR1 (ICOUNT) 

IUR=LP?(ICOUNT) 

IKI.fO.O)  370,300 
IF  ( I . EO.  I'JQ  + 1 )  290,310 

F=(PMIN-(2.0»DELF) ) 

J=J*1 

AKi=(A.2fl31853*r)/CT 
AK2=(6,?531ft?3<>F)/C2M7N 
FO= !39.47eAl76»F*F) 

IGV=I7FRO=0 

CALL  ITERATE (AK) , FD , H ,C2D,DLS , HD,HD3 , CTD , I T ,NCR*R , Z, A , i t ,DL , I ZERC 
»l  > 

IFl’NCO.GE.I)  330,340 

CALL  HALF  (AK1,AK2,FD,NCR,A,I.,  IC) 

DKN(IC,J)=AKN(IC) 

C(IC,J)= (6.2831 «53*F)/DKN(IC»J) 

GO  TO  350 
C'IC«J)=0.0 

f=f+oflf 

IMJ.EO.K)  360,320 
CALL  0«QFRROP(0,0> 

1  =  1*1  & 

IC=IC*1 

J=0. 

GO  TO  300 

L=*-2 

ICOUNT=0 

IC=1 

J=2 

RFI.AG=O,0 
ICOUNT=ICOUNT+1 
■I=L<31  (ICOUNT) 

IUP=LR?(ICOUNT> 

IF(I.fo.O)  520,390 
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>5 


10 


n 
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390  J=J+1 

IF(C<IC.J-?) .En.0.0)  400*450 

400  IF(C(IC.J-1) .EO.O.O)  410.430 
410  IFJCdC.J)  .EO.O.O)  4?0.440 
280  420  GV ( IC. J-25 =10.0510 

GO  TO  470 

430  CPUC,J-2>  =  (CdC»J*l)-CdC.J-l))/(12.5663707*DF|_F> 

GO  TO  460 

440  CPdC.J-2)  =  (-CdC.J«-2)*4.0«C<IC»J*l)-3.0»C(IC.J>  >/(12.5663707*DEl 
*) 


60  TO  460 

450  CPdC.J-2)  =  (-CdC«J  +  ?>+n.0»CdC.Jd>-A.0«CdC.J-l>+CnC.J-2))/ 
«■  <75.3962235*DELF) 

460  VdC*J-2)  =  d.0/CdC»J)  >-(OKM(IC.j;/C(lC. J)  )*CP(IC.J-2) 
GV(IC.J-2)  =  {1.0/v;iC.J-2) ) 

470  IF(J.FO.L)  400.390 
460  MMIJ  (10=1 
IC=IC+1 
1  =  1  +  1 
J=2 

IF ( IC.FO, 13)  490.510 

490  PRINT  16000 
IK=K~1 
PRINT  25000 
300  PRI NT  11000 

PRIMT  20000.  (MM(j(J)  ,J=1,  IK) 

PRINT  26000. LFRF, (LGPU.  J2=1.IK5 
F=FMIN 


DO  500  N=3,L 

PRINT  22000, F,  (GV(J,M-2.',J=1,IK) 

500  F=F+nFLF 

IF (°FLAG.FO,2.0)  GO  TO  530 
510  IFd.ro. IU9+1  )  360,390 

520  RFL  46”? . 0 
GO  TO  490 

C  SU^ROtJriNF  PLOTS  GROUP  VELOCITIES 
530  IF  (GRIIPLT  .  fc’Q. ] . 0 )  540,9997 

540  CALL  GPOijP  (  IC.K.FMIN.FMAX.DELF) 
9997  CONTINUE 

415  999 6  CALL  PLOT ( 20 . , 0 . , 3 » 1 ) 

9999  CONTINUE 
END 
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30 


35 


40 


45 


30 


35 
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30 
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SUBROUTINE  GROUP(IC.K,FMIN,FMAX,DELF) 

COMMON/ A A/CP (20,20) , V ( 20 .20 ) ,GV (20 .20) ,C (20  *  20) »DKN (20*20) 
KK-K-4 

CALCULATION  TO  FINO  “AX.  AMO  MIN.  GROUP  VELOCITIES 
GV.MAX  =  0.0 
00  20  1  =  1  * IC 
DO  10  J= 1 . KK 

10  I F  ( G V  ( I ,  J )  »GE  .GVMAX.  A.NO.GV  ( I  »  J)  ,N£  .  1  0. 0*«  10)  GVMAX=GV  ( I  ,  J) 
20  CONTINUE 
GVMIN=GVmAX 
DO  40  1=1. IC 
00  30  J= 1 . KK 

30  IF (GV(I.J) .LE.GVMlN)  GVMIN=GV (I  * J) 

40  CONTINUE 

SCALING  OF  GROUP  VELOCITIES 


IGVMIN=GVMlN/10.0 

GVMJN=IGVMIN*10.0 

GVOIFF=GVMAX-GVmIN 

IF (GVD I FF.LE. 10.0) 

50.60 

SO 

6Vm 1=1.0 

SCALE=1.0 

GO  TO  120 

60 

IF (GVD IFF. LE. 20,0) 

70.80 

70 

G VN I =2 . 0 

SCALF=2.0 

GO  TO  120 

1 

IF (GVOIFF.LE.40.0) 

90.100 

90 

GVM 1=4,0 

SCALE=4.0 

GO  TO  120 

100 

IF (GVOIFF.LE.60,0) 

.0,115 

110 

GVM I =6.0 

SCALE-o.O 

GO  to  ]20 

ns 

IF (GVD IFF .  LE.80, 0) 

1 16.120 

110 

GV*-'  1=0.0 

SC«LF=P.O 

SCALING  OF  FREQUENCIES 


120 

FOIFF=FNAX-FMIN 

IF  (FOIFF.LE.IO.O) 

130, 

140 

130 

FNI=1 .0 

FSCALE=1 .0 

GO  TO  260 

1  4  0 

IF (FOIFF.LE.50.0) 

150, 

160 

150 

FNTsS.O 

FSCALE=S.O 

GO  TO  260 

160 

IF(FOIFF.LE. 100.0) 

170 

O 

CO 

► 

170 

FN 1=10,0 

FSCALF=] 0.0 

GO  TO  260 

IF  (FOIFF.LE. 250.0) 

190 

,200 

190 

FN I =25 . 0 
FSCALF=25.0 
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65 
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r5 


^5 
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5 
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GO  TO  260 

200 

IF  (FOIFF.LE. 500.0) 

210. 

220 

210 

FNI=*0.0 

FSCALF=50.0 

GO  TO  260 

220 

IF  (FOIFF.LE. 1000. 01 

230 

,240 

230 

FNI=lO0.0 

FSCALE=100.0 

GO  TO  260 

240 

IF(F0IFF.LE.15C0.0) 

250 

,260 

250 

fnt=igo.o 

FSCALE= 1 60 , 0 

C  PLOT  PACKAGE  FOR  GROU°  VELOCITIES 


?6o  call  plot<s.5.-io.,3.i) 

CALL  AXIS (0.0, 0.0, 21HGPOUP  VELOCI TY < M/SEC ) *  21 , 1 0 . 0 , 90 . 0 , GVMIN, 
*GVMI) 

CALL  AXIStO. 0,0.0, 13 HF RE OUENCY(HZ) ,-13,10. 0*0.0,  FMIN.FNI) 
00  350  1  =  1  ^IC 
OUt-'F =FMI  N 
00  340  J= 1 i KK 
X= (OUMF-FMIN) /FSCALE 
IE(GV(J,J).FO. 10. 0**10)  260  *  290 

?«0  nuwF  =  nii^F  +  OELF 

CALL  PLOT { X , Y , 3  »-l ) 

30  GO  TO  3^0 

2P0  Y=(GV(I,J)-GVMIN)/SCALE 

I F (GV ( I. J-1).F0.1 0.0**101  294.295 

29a  CALL  PLOT (X.Y. 0,-1 ) 

2Q6  IFCJ.PQ.l)  300.310 
300  CALL  PLOT (X.Y, 3,-1) 

GO  TO  320 

310  CALL  “LOT (X.Y, 2,-1 ) 

320  OUMF  =  DMMF  +  OELF 

IF (J.FO.KK)  330,340 
330  X=X»0.1 
Y=Y-0 . 03 

CALL  SYM9L4 ( X , Y , 0 . 07 , ahmODE , 0 . 0 , 4 ) 

X  =  X  +  0.6 

xim=float  <  i  > 

CALL  NllMMP(X,Y,0.O7,XlM,0.n,-l) 

340  COK'TIMUF 
350  CONTINUE 

CALL  PLOTU5. 0,0. 0,3,1) 

360  RETURN 
100  END' 
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SUBROUTINE  HALF  COC  6400  FTM  V3.0-P316  0PT=1 

SURPOUTINE  HALF ( AK) . AK2.FD.NCR. A , I , IC) 

COmmon/1/  XT (  3100) 

COmmon/2/Z2N ( 3 1 00 ) ♦  C  2  0 ( 3 1 0  0 ) 

5  5  C0MMOH/Afl/CP(20,20)  t V (20.20 ) .GV (20.20 ) .C (20.20 > .DKN (20.20) 

CO^MOM/AR/Zl  (SO)  » C  1  (GO )  *  Z 1 N ( GO )  .LLB(IO)  .LUB (10)  .  MflU  ( 30  )  .  LP1  (10)  . 
*L«? (10) 

COMMOM/AC/AKN (  ]2).RR(  1?>.PV(  1?),AM(  12).ZF(12) 

COMMON/ An/EPS ILON,h,h0.ho3.0L.DLS.R. I T.CT.CTD.MODSHAP  .MODEPLT. 
10  *6«(JVR|.,fiRUPLT.  I7FP0*  JGV 

EQUIVALENCE ( XT « ZZ ) 

1000  FOPmat ( 0  UPPER  AMPLITUDES  ARE  7ER0ED  FOR  THIS  MODE  *> 

C  HALF  INTERVAL  SEARCH  FOR  THE  knS  ASSOCIATED  WITH  THE  MOOES  DESIRED 
JUMP=0 

IS  10  AKI=AK?-AK1 

OIVIOF=2.0 

11  BKN=AK) ♦ (AKI/DIVIOE) 

IZFR0=0 

CALL  ITERATE (8<N,FD.H.C2D.DLS.HD.HD3»CTD» IT.NCR.R.Z. A.ZZ.OL. IZEP 


20  *1) 


20 

IF(JUMP.EO.l)  40.20 
IF(MCR.FO.I-l)  50.30 

30 

IF(NCo.LT.I-l)  35.36 

35 

DIVIDF=niVIOE+1.0 

RE 

GO  10  11 

36 

AK 1 =  GKN 

40 

GO  TO  10 

IF (MCR.EO. I )  70.A5 

45 

IF(NCR.LT.I)  35.36 

30 

30  50 

JlJMPrJUMP.l 

60 

AK?=HKN 

70 

ZT 1  =Z 

GO  TO  10 

IF  (Z)  80.160.90 

35 

PO 

aknl=«kn 

iO 

90 

ZL  =  Z 

AKNP=AK2 

ZR=7T1 

GO  TO  98 

AKN»=RKN 

ZR  =  Z 

AKNL=AK2 

ZL=ZT 1 

98  CONTINUE 

•  5  C - : - 

C---A00EO  3/?  HALVING  APPROACH 

IF  (ARS(ZL) .LT.l .ES0.OR.ZR.LT.1.E50)  GO  TO  100 
DIFF?= ( AKNR-AKNL) /?. 

AKA=AKNL+0IFF2 
.0  101  TZFRn=o 

CALL  ITERATE (AKA.FD.H.C2D.DLS.HD.HD3.CTD. IT .NCR. R. Z . A . ZZ.DL ♦ IZERf 

«.I) 

IF  ( ARS  ( 7 ) .LE.EDSILON)  170.9) 

91  IF (7P*7) 92.94.94 
'S  55  9?  DIFF=(AKNR-AKA)/2. 
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60 


6? 


70 


75 


M  O' 


*0 


IFIARSIDIFF) .LE.l.E-1?)  160,93 
93  AKNl=AKA 
ZL=Z 

aka=aka*diff 

60  GO  TO  101 

04  OIFF= (AKA-AKNL) /?. 

IF(APStOIFF)  .LE.l.E-1?)  160,95 
95  AKNR=AKA 

AKA=AKA-OIFF 
ZP  =  Z 

GO  TO  101 

0— - 

C  LOCATING  THF  KN  ASSOCIATFn  WITH  THE  MODE  DESIRED  USING  THE  KNL  AND  K 
C  FOUND  AROVE  RY  THE  METHOD  OF  FALSE  POSITION 

100  AKA-AKNL+ (ZL° (AKNR-AKNl) ) / ( ZL-ZR ) 

1001  FORMAT(I5,RF\6.10) 

IZFR0=0 

CALL  ITERATE  ! AK A , FD , H, C?D» DLS » HD , HD3 , CTD, I T ,NCR . R* Z, A, ZZ , OL *  I ZER 
*1) 

IF  ( A P. S  (Z)  .LE- EPSILON)  170,110 
110  IF ( ?P«? )  120,140,140 

120  IF ( A R S  (AKNL-AKA) .LF.1.0E-12)  160,130 

130  ZL=7 

AKNL=AKA 
80  GO  TO  100 

140  IF ( A  H  S  (AKNR-AKA) .LE.l .OE-12)  160,160 

150  ZR=Z 

AKNP=AKA 
GO  TO  100 
160  IZ^PCcl 

PRINT  1000 

165  CALL  ITERATE (AKA,rn,H,C?D,DLS,HD,HD3,CTD, IT,NCR,R,Z, A,ZZ,DL, IZERC 
«I) 

170  AKN(IC)=AKA 
ZF(TC)=Z 
I  GV  =  0 
RETURN 
100  ENO 
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CDC  6400  FTN  V3.0-P316  OPT=l 


SUBROUTINE  IT£RATE(AKN,FD,H.C2Q*DLS,hD,HD3*CTD,IT,NCR*R,Z* A,ZZ*' 
«IZFRO, I ) 

dimension  zztiT)  ,c?D(it> 

5  5  C  SUBROUTINE  USES  A  FINITE  DIFFERENCE  TECHNIQUE  TO  NUMERICALLY  FIND  Tl 

C  SHAPE  ASSOCIATED  WITH  A  GIVEN  KN 

C  CONSTANT  IN  FRONT  OF  FUNCTION  ASSUMED  TO  BE  ONE  INITIALLY 
MAXFLAG=0 
N=IT 

10  A=KD» ARS  (AKN«AKN-(FD/CTD) ) 

ZM1=R 

XND=HD° ( (FD/C20 (IT)) - ( AKN^AKN) ) 

XNU=HD* ( (FD/C20 (IT-3 ) ) -<AKN*4KN) ) 

C - next  CARD  ADDED  AS  TEST  NOL  3/6/73 

1*>  XNW=Hn*(  (  (FO«  (SORT  (C20  ( IT-1 )  )  )  > / (C2D ( I T> *»1 .5)  )  -  ( AKN^AKN)  ) 

C  TAYLOR  SERIES  EXPANSION  TO  OBTAIN  NEXT  Z 

Z=ZM1* ( 1 . + ( ( DL/R ) *SOPT (A) ) -(5.»DLS»XND/6.) + (DLS*XNW/3. ) - 
* ( ( (DL«DLS)/6.0) *SORT (A) *XND* { 1 ,/R »  > ) 

NCR  =  0 

20  IF(ZMHZ.LE.O.O)  10,20 

10  NCR=NCR+1 


20  PS=(?.0-(HO«OLS)<U  (FD/C2D ( IT-1 ) )~(AKN*AKN) )) 
ZP1=(PS*Z)-ZM1 
IF  (  Z  *  Z  P 1  .LE.O.-O)  30,40 


2' 

30 

30  NCR=NCR* 1 

40  ZZ ( I T ) =ZM1 

ZM1=Z 

Z  =  ZP1 

C****CHANGE  NOL  4/6  REMOVE  DOWN  TO  77 

30 

AO  I T  =  I T  —  1 

IF ( IT-2, LT  .0)  90,20 

90  I T=N 

ZZ(?)=ZM1 

ZZ(1) =7 

35 

GO  TO  130 

100  I  T= I T- 1 

IF ( I T-2 . LT . 0 )  120,110 

110  ZZ ( I T ) =0 . 0 
,  GO  TO  100 

40 

120  IT=N 

ZZ  ( 1 ) =7=0.0 

ZZ  (?) =7M1  =  0,0 

130  COMTINUF 

C**«mh>chanCF  NOL  4/6  ADD  TO  140 

45 

IF ( I  ZERO. N£ . 1 )  GO  TO  140 

IF (NCP»GT  » { 1-1 ) )  131,136 

131  DO  134  IT=1,N 

IF(  (77(IT)<*ZZ(ITM) )  .LF.O.)  GO  TO  135 
Z7 ( IT) =0. 

5C 

134  CONTINUE 

135  Z7 (IT) =0. 

GO  TO  1139 

136  DO  139  IT=1,N 

IF ( ABS(ZZ ( 1T+1 ) ) ,GT. ARStZZ ( IT) ) )  GO  TO 

55 

55 

ZZ ( I T ) =0 . 

NOLTR  7U-95 


60 


IWOQUTINF  ITERATE  CO'C  6400  H'TM  V3.0-P316  0PT  =  1 

139  CONTINUE 
1139  IT=N 

140  CONTINUE 
RETURN 

60  ENO 
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CDC  6400  FTN  V3.0-P316  OPT  =  l 


5 


10 


IS 


20 


30 


35 


40 


45 


SUBROUTINE  SHAPE(N.I,IMODE) 

CO^MON/1/  XT (  3100) 

CObmON/?/Z2N(3100) .  C2D  (3100) 

5  COMMON/ A  4/CP (20.20)  ♦  V ( 20 , 20  >  *GV  (20*20)  ,C (20 , 20 )  .D*N(20»20) 

COMMON/ AB/7.1  (50)  .Cl  (50)  ,Z1N(50)  ,LLP(10)  ,LU3(10)  ,MMU(30)  .LP1  (10)  , 
»l«?(10> 

C0MM0N/AC/4KN (  1?).PR{  1?).PV(  12)  t  AN (  12) e  ZF ( 12) 

COMMON/ AO/EPS  I  LON, h. HO. HD3. OL.DLS.R, I T.CT. CTO. MOOSHAP  .MOOEPLT » 
»GRlJVFL.G»UPLT,  I7t'RO,ir,v 
Common/ A E/LDEPTh  .LAMP  . LGRU  .LFRE 

DIMENSION  ZT (3100) .77(3100) 

EQUIVALENCE (22N.2T) . (XT.ZZ) 

C  I  PERTAINS  TO  IC  IN  MAIN  PROGRAM 
DO  in  JT  =  1 » N 
XT (  JT)=ZZ(JT)»AN(I) 

10  ZT ( JT) =72N< JT) 

WRITE (R. 12) I  MODE 

WRITF (R, 14) (ZT (LL) .XT(LL) ,LL=1,N) 

12  FORMAT (110) 

14  FORMAT (4E20.10) 

IF(MOnFPLT.EO.l)  20.40 
C  PLOT  PACKAGE  FOP  MOOE  SHAPES 
20  I F ( I ?FPO»EG. 1 )  24,25 

24  XB=6.5 
XD=-2.0 
XL=4.0 
XI=-4.0 


GO  TO  26 
30  25  XG=4.5 

XO=-1.0 
XL=2. 0 
X  I  =  -2 . 0 

26  CALL  DLOT (  XB,0.0  .3.1) 

CALL  A X I S  (  XD , n . 0  « 20HNORM AL  I  ZFD  AMPL  I  TUOE  .  20  .  >'L  ♦  0 , 0  ,  XI  ♦?.  0  ) 
CALL  AXIS(0. 0,-6. 0.3  4HNORM  AL  I  ZED  DEPTH.  1  .< ,  5 . 0  ,  QO  .  0  ,  <  l»0»-0.2) 

CALL  SYMAL4 (-0.6,-6.0,0.096, 4HM00E ,0.0,4) 

XIR  =  FLOAT ( I  MODE ) 

CALL  NUMSR(0.26.-6.0,0.095,XIB,n.O,-l ) 

NN=3 

00  30  JT=1,N 
X  =  X  T  (  JT)/2.0 

Y=-ZT ( JT) 05.0 
CALL  PLOT ( X , Y , NN , - 1 ) 

^5  30  NN=? 

40  RETURN 
END 
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COC  6400  FTN  V3.0-P316  OPT-1 


5 


10 


15 


30 


35 


AO 


45 


5Q 


5$ 


5 


c 


c 


c 


SUBROUTINE  ZCINTpc  (NP«N»E »C2MIN»PR0pLT) 

CO^mOm/I/  XT (  3100) 

COmmON/2/Z2N(3100) .  C200100) 

COMMON/ A  A/CP <?0, 20)  ,V  (20*20)  .GV  (20,20)  ,C  (20,20)  ,0^(20,20) 
COMMON/ am/ZI  (50)  .Cl  (50)  *7 IN (50)  ,LL8(10)  »LU8(10)  ,MMJ(30)  ,LH1 (10) 
«L«2(10) 

COMMOM/AC/AKN  (  12).PP(  12)  ♦  P  V  (  1?)»A\'(  1?)  *ZF  (12) 
C0MM0M/An/5PSILON,H,HD,HD3,0L*0L5f Pt IT*CT»CTO«MQOSHAP  » MODEPLT , 
*GRUVFL  , GRU°LT  <  I  /ERG, IGV 
COMmom/ae/LD£PTh  .LAMP  tLGRU  *LFPE 
LINPAR  INTERPOLATION 
DO  50  J= 1 « N 
00  40  1=1. NO 

IF  ( 7 1 N  ( I )  .EO.Z2N'(J)  )  10,20 
10  C?D(J)=C1(I) 

GO  TO  50 

20  IF(ZlH(I).GT.Z2MiJ))  30.40 
30  ZGT=Z1M(I) 

ZLT=Z 1 N ( I -1 ) 

CGT=C1  (I) 

CLT=C 1 ( I -1 ) 

EQUATION  USED  FOR  LINEAR  INTERPOLATING 

C20(J)=(<Z2N(J)-ZLT)/(ZGT-ZLT))*(CGT-CLT)  «•  CLT 
GO  TO  50 
40  CONTINUE 
50  CONTINUE 

SEARCH  FOR  minimum  ROUND  VELOCITY  ON  SOUND  VELOCITY  PROFILE  (ALSO  F 
MAXIMUM  SOUND  VELOCITY) 

C2ma*=0.0 

DO  AO  1=1. N 

60  IF (C?0 ( I > .GE.C2MAX)  C2UAX=C2D(I> 

C2MIN=C2mAX 
DO  70  1=1, N 

70  IF (C20 ( I ) .LE.C2MIA )  r?MlN=C2DCI) 

PLOT  PACKAGE  FOR  SOUND  SPFED  PROFILE 
IF (PRnPLT .EO. 1.0)  75  » c  C  0 

75  CDJFF=C?maX-C?min 
I C2M=C?M IN/10,0 
C2M= I C2M&  10,0 
IF (CDIFF.LE.5.0)  80,90 
80  PL =2.5 
PI=0.5 
C2m 1=1.0 
SCa'lF  =  0.5 
GO  TO  170 

90  IF  (C f' IFF, LE. 10.0)  100,110 

100  RL  =  2 . 5 
P I =0 . 5 
C?MIr?.0 
SCALE=0.25 
GO  to  170 

110  IF  trDIFF,LE.25.0)  120,130 
120  pL=2.5 
PI =0 . 5 
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SUBROUTINE  7CIMTE« 


COC  6400  FTN  V3.0-P316  OPT=l 


60 


6e* 


70 


75 


Rr 


°5 


60 


C2”I=5.0 
SCALE=0.1 
60  TO  170 

130  IF (C0IFF.LE.50.0)  140,150 

140  PL=5.0 
PI  =  0.5 
C2”I=S.O 
SC4Lf=0.1 
60  TO  170 

150  IF<CDJFF.L£. 100.0)  160,170 

160  PL=5.0 
PI  =  0.5 
C2“I=10.0 
SCALE =0 . 05 

170  CALL  PLOT (0. ,0. *3, 1 ) 

PL=PL+PI 

CALL  AXIS(0.0,-S.0,16HNORM4LIZED  DEPTH.  16,5.0,00.0, 
CALL  AyiS(0.0,0.0,lAwSOUND  SPEF.O  (M/SEC)  ,  In, PL  ,0.0, 
NN=Z 

00  IPO  1=1, N 
X=(C?D(I)-C?M)#SCALE 
Y=-72N ( I ) *5 . 0 
CALL  PL0T(X,Y,NN,-1) 

1R0  NN=? 

CALL  PLOT (2. 5,0. »3,1 ) 

200  CO.NTINIIF 

DO  210  J=1,N 

210  C?D(J)=C20(J)*C2D(J) 

RETURN 

END 


1.0. -0.2) 

C2M , C2M 1 5 
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SIMPOUTINF  LINK 


CDC  6400  FTM  V3.0-P316  OPT=l 


SUPPOUTINK  LINK  (X»Y*N,K) 
OIMKNS I  ON  X(l) ♦  Y(l) 

13=3 

NP=N«K 

5  oo  in  i=ifNp,K 

CALL  PLOT  (X(I) *Y(I) »I3»-1> 
10  13=? 

UPTURN 

END 
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Appendix  lilt  Normal  Mode  Summing  Program 

PROGRAM  NORMSUM  COC  6400  FTN  V3.0-P 


B 


10 


15 


20 


25 


30 


35 


40 


45 


50 


PROGRAM  NORMSUM ( INPUT, OUTPUT, TAPE7  »TAPE1) 


TAPE7  IS  OUTPUT  PILE  OR  TAPE  FROM  NORMAL  MODE  PROG. 

SR-  STARTING  RANGE  FOR  CALC. 

RINC-RANGE  INCREMENT 
NR-NUMPER  OF  INCREMENTS 

PRO (1.2.. ...5)  ARE  UP  TO  5  RECEIVER  DEPTHS 
SOURCE  IS  SOURCE  DEPTH 

ABOVE  MUST  HE  ONE  OF  DEPTHS  USED  IN  NORMAL  MODE  PROG. 

NPRTN  .GT.  0  —  PRINT  RESULTS 
NPLOT  .GT.  0.  — °LOT  RESULTS 
PRINTR-INITI AL  RANGE  ON  PLOT 
PRTNC-RANGE  INCPEMENT/INCH  ON  PLOT 
PODINT- INITIAL  PROP  LOSS  AT  ROTTOM  OP  PLOT 
PD9IN-  PROP  LOSS  IMCRFMENT/IMCH 

CARO  RGAiJ  •STATEAlEKJTS- uiNgs  3  e,  3  1 ,  33  ,36" 

COMPLEX  TEST 

COMPLEX  SUM1.PS.PS.  Z.CSI4 

OlMENSIONAKNf  1000) .  XTS(IOOO) i XTR ( 1 000  • 

DIMENSION  PRD  ( 5 )  ,  TS  U  000 )  .  FS  (500 )  »  Tl.  <  500  )  .  RO  (S )  »  AXT  ( 1 2 ) 
DIMENSION  IP(S) .TPFNO(S) 

CALL  PLOTS(0.,0.,1) 

READ  1000.SPiRIMC.NP 
1000  FORMAT (?Flo. 3. IS) 

30  IF (NR.GT.500)  STOP  6 

RF AO  1  020,  ( PRO ( I L ) . IL=1 »5) 

1020  F  ORM at ( RF 1 0 . 3 ) 

READ  1021, SOURCE 

1021  FOPMAT (F 1 0 . 3 ) 

READ  1010 ,NPPIN,NPL0T.PRINTR.PRINC,PDRINT,PDPTN 
1010  format (215, 4F in. B) 

RE AD (7, 7 600) SOEP, <RO( IL) , IL=1 ,S) ,LI 
7600  FORMAT ( 6F 1 0 . 6 , IE) 

WETEST=. 5/FLOAT (LI) 

REAO(7,11SOO)F. M 1 , H, MCP 
NCP1 =NCR 

11500  FOPMAT (3F10.3,I4) 

DO  7  I=1,S 
IR ( I ) =1 

7  PRO ( I ) =PRO ( I ) /H 
SOURCF=SOURCE/m 
PJ=3. 14159265359 
CS 1 4  =  CMPLX (0. ,-p;/A. , 

SUU 1 = (p ] /h) °SQRT (2.0«3. 1A1592)°CFXP(CSI4) 

J=0 
IS=1 

10  J=J*1 

READ (7, 2000) AKM( J) ,PV,RR 
2000  FORMAT (3E20. 13) 
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P^'-RAM  NORMSUM  COC  6400  FTN  V3.0-P316  OPT-1  C 

IF (EOF (7) ) 20, 10 


20 

J=J-1 

1 

25 

READ (7 *3000) AOFP. (A  XT (III) ,111  =  1,12) ♦IFl 

> 

3000 

FORMAT (F10.6, 3^20. 13/4E20.13/4E20.13/E20. 13* 12) 

60 

IF(E0F(7) 130.36 

30 

J=J+1 

READ (7,2000) AKN(J) *PV«RR 

IF (EOF (7) ) 60, 10 

35 

IF(IFl.EQ.O)  GO  TO  45 

5  00  37  IJK=1,S 

IF  <  (IF1 .ER.l) .AND. (ARS(RO(IJK)-AOEP) .LT.WETEST) .AMO. (RD ( I JK) .NE. 0 
*) )  60  TO  36 
37  CONTINUE 
STOP  6 

0  3fi  IRFNO ( I JK) =IW ( I JK) +1 1 

ITP=IPENO(IJK) 

ITR=IP (IJK) 

ICO=0 

00  3°  IP=ITR. ITP 

6  IC0=IC0+1 

39  XTR(IP,IJK)=AXT(ICO) 

IR(TJK)=IR(IJK) +1? 

GO  TO  ?6 

45  IEMD=!S+ 1 1 
10  ICO=0 

DO  47  I JS=IS« IE-NQ 
IC0=IC0+1 

47  XTS ( I JS ) = AXT (ICO) 

IS=I 6* 1 2 

*6  85  GO  TO  25 

60  IF(A5S(S0UPCE-S0EP) .LT.WETEST)  GO  TO  64 
OG  6?  IC11=1,5 

IF (ARS ( SOURCE -R0( IC11)). LT.WETEST)  GO  TO  63 

62  CONTINUE 

RO  STOP  3 

63  00  1061  IJJl=1.100O 


1061 

XTS(IJJl) =XTR(IJJ1 , IC1 1 ) 

CONTINUE 

i 

S 

t 

64 

00  000  I C=  1  * 5 

r 

96 

IF(PPO(IC) .EO.O.)  STOP  ? 

5 

1 

65 

R=SR 

DO  66  TO=l , 10O0, 

TS(lO)=0. 

i 

i 

00 

66 

00  66  101=1, 5A0 

FSU01)=CMPLX  (0.,0d 

1 

3 

| 

OO 

00  PO  I J= 1 , 6 

IF ! ARS (DR0 ( IC ) -P0 ( I J) ) .LT. wETFST)  GO  TO  96 

CONTINUE 

\ 

C - 1 J 

IS  INDICATOR  FOR  XTR;1000*IJ>  ARRAY 

(. 

OR 

00  inn  im=i,mcni 

100 

TS(  JM)  =XTS(  I--*)  «XTP(  IN.  I  j)  /SORT  (  AKN(  IN)  ) 

1 

110 

110 

no  200  jj=i,mr 

PS=  Cpplx (0.  ,n. ) 

DO  1«0  I J J= 1 , NCR  1 

Z=Akn(ijj><»R«  Cnplx  (n, ,  1 , ) 
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CDC  6400  FTN  V3.0-P316  OPT=l  0 


180  PS=PS+TS(IJJ)*CEXP(7) 

FS(JJ)=FS(JJ)+(SUM1 /SORT (R) )  *PS 
TL (JJ)=-20.*ALOG10 (CAHS(FS(JJ) ) ) 


t 


0 


5 


0 


5 


.0 


55 


115  TL  ( JJ)=TL  <JJ)  ♦?0-*ALOr-10  <  1./.015A) 

C - CONVERTS  FROM  RE  1  PETER  TO  RE  1  YARD 

200  R=R+RIMC 

RA=PRD(IC)»H 
SA=SnURCE»H 
PRINT  2500  .  S  A  .  P  A 

IF (NPLOT.EO.O  )  GO  TO  400 
CALL  PLOT (2. .0. ,3. 1 ) 

CALL  SYM8L4 (0.. 3 .*.14.1 3HS0URCE  DEPTH=,90. ♦ 13) 

CALL  NUMBR(0. .5.2.. I4.fi, 90. .1) 

CALL  PLOT ( 1 • »  0 • »  3 • 1 ) 

CALL  SYPBL4(0. 14. 15H»eCEIVER  DEPTri=,90. . 15) 

CALL  NUMFjR(0.,5.5..14,RA,90.»1) 

CALL  PLOT (1..0..3.1) 

CALL  SYMRL4 (0. , 3. • .14. 1 4HFREQUENCY (HZ ) = .90 . ♦ 1 4 ) 

CALL  NUM5R(0.i5.5».14,F.90..“l) 

CALL  PLOT (2..0..3.1) 

CALL  AXIS ( 0. , 0. .RPR A NOE (KM) .-9, 10. .0 . .PRINTR.PRINC) 

CALL  AXIS(O..Q..20PPROPAGATION  LOSS (OB) . 20 . 9. » 90 . .PDRI NT. PDRIN) 
CALL  PLOT (0. ,0. .3.-1 ! 

400  CONTI  MUF 

2500  FORMAT (?8X.oN0°PAL  MODE  PROPAGATION  LOSS  RE  1  YD*//25X » ^SOURCE  DE 
*TH=»,F10.3.10X.*RECFIVER  DEPTH  =  «  .  F 1  0 . 3//30X .  <»  RANGE  *,«TRANSMI 
140  *SION  LOSS ( DtO  0 ) 

R  =  SR 
NN=3 

^0  600  I MPs 1 . NR 
IF (MPRIM.EO.O)  GO  TO  GOO 
PRINT  3500,R,TL(TPP' 

500  IF(NPLOT.K.Q.O)  GO  TO  600 

X=((R/(1.E  +  3)1 -PRTNTR) /PR  INC 
Y=  (POR1NT-TL ( I mp ) ) /APS (PDBIN) 

Y=AMAX1 (0. ,Y) 

CALL  PLOT (X,Y,NN.-1) 

NN=? 

600  R  =  R*R  T  NC 

3500  FORMAT (25X.F10, 3.  5X.F10.3) 

CALL  PLOT (li.,0.,G.l) 

155  800  CONTINUE 

END 
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Appendix  tv,  Evaluation  of  Second  Range  Derivative 


In  order  to  evaluate  pressure  on  the  caustic  (Equation  4.1), 
we  need  to  evaluate  &  u/DC.1,  and  obtain  S  (Equation  4.15). 

is  defined  in  terms  of  the  derivative  of  an  integral t 


*  5  L 


(IV. 1) 


But  the  integral  can  be  split  up  into  several  integrals,  each 
evaluated  in  a  different  sound  velocity  profile  layer: 


l 

3Wc  - 


. .  r  f4**'  a*l 

1  \  \  (iv.2) 

ifc  a*  i{  ^§’5* 


where  ma)  * 

Thus  we  can  evaluate  or  ,  in  each  layer  and  then 


sum  to  obtain 

In  CONGRATS,  A  R  in  a  layer  from  to  z2  is  defined  as  (39): 


R  :  l4Rl  * 

fo  +  s^*)  + 

Gy  / 

Pft(ro-1)]? 

(iv. 3) 

where 

J 

42  r 

42,  ♦  62L  i 

■& 

42: :  i  -  2-ft 

(IV.3A) 

D  5 

2»-24  > 

(IV. 3B) 

X«  :  C  A2ft  +b 


v 


eu-b4 
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-1 


a  :  V,  -  Cy 

k  :  aSt  ‘43. 

C. 

Is  =  6 -5»'>  8  3,  ‘  4  3.3* 


H  is  one  of  several  results  depending  on  the  sign  of  c  and  q. 
g q,  g1(  g2  and  vQ  are  the  four  parameters  (Equation  3.1  ) 
in  the  particular  layer  being  considered. 

Weinberg  then  obtains  Sty  DC*  ,  the  derivative  necessary 
for  amplitude  calculation  along  a  rayi 


.  / 1  ^0  .1  \r  ♦  zt>  f  q  'Daa  ♦  b  fi/ii  -jm 

5cv  -  (  d  cv  r  r  13;  k  p  jJj(iv-4> 


He  then  evaluates  the  necessary  derivatives! 


!>2  .  ,  and  2 

la.  »  '  a  '  *s  i 


*dcv 


"ac 


[<?.  -■>* 


(IV. 5) 


for  the  possible  vaiues  of  c  and  q,  i.e.  positive,  negative,  or 
aero.  For  calculation  on  a  caustic,  we  need  one  moie  derivative, 
V’R/t )Cy  •  Taking  the  derivative  of  Equation  (IV. 4)  yields i 

1 R  ♦  £  ♦  JL  *1  00  *&R 

'DC*1  D1L^CyJ  D  7)cvl  cw%  D  0Cy  -SCv/ 


-J.  ’DR 
Cv  ^C14 


*  f  i.  ?  f  3  ("j.  /  H  _  |\]} 

tcvxv  c7*j[^ac,  prc,[t la'll 


*aoq,  ■££* 

TJ 


■sop  t  ri/H 

I  ^TUo  'J 


(IV. 6) 
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NOLTR  7^-9^ 


New  derivatives  required  are: 

~£j&  ^ 

7)  C„x 


SL  fi  (h  -  ,)' 


(IV. 7) 


Bt 

For  c /  0 


*£& 

7>CVV 


-  C> 


$&  * 

acv  2[v,/w  vt,/l]  1  V4 


f  ^Y</2Cy  +  ’dVac 

V1 


(IV. 8) 


Then 


•£d 

*t>C/ 


_g _  f  ^./3c.  .  'aVac, 

[.  "v/» 


-  D 


f  AaoVj^-i 

1>  y*  2Y,^ 


*[V*‘V*] 

So  we  also  need  ^ 

For  c  =  0  (Ray  vertexes): 

^  r  O  ,  Vt  SO 

Therefore  terms  multiplied  by  •  ,  or  by 

i  —  ^Cv 

3  a  a 

Cl  tiCg1 

From  Equation  (IV.3A) 


Vf 

?c„ 


(IV. 9) 

rxf? 


l>Cg 


t 

z  2CV 


3jjT4  r  j.  £ 

1  *cv 


(IV. 10) 


are  zero. 


(IV. 11) 
(IV. 12) 
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The  quantity  ‘'Li,1;  is  non-zero  on3y  in  the  layer  in  which 
'a<*t 

the  ray  turns.  Only  here  will  z2  be  a  function  of  cy. 

D  ■■>•••  _ 

1  'iCy'- 

Prom  Equation  (IV.3B) 

f. s 

except  for  turning  point  layer  where 

-  0  _  (J1 

Then  for  a  non  vertex  layer  (  ~r*  £  °) 

^  y  Cy  / 

=5c>  ^ 

For  a  vertex  layer 

eJjU  -  ~  (I  +  a  Ait)1 

«v  «  *<  -  —  n  /  T1 


Since 


+  ( I  +  5,  2f?  . 

v  *acv 


3a  ;  Z  C,' 

•SCV 


'DCv1 


i  Fi  1  h  -j)l 

El  t>Cvl  [dlD  'J 

Starting  with  £  J.  ^  ^  _  , )  “J  ,  taking  tv.-o  derivative 
recombining  several  terms,  it  can  be  shown  that 

♦ 

fi  ( «  „i)l  -  -2.  [2  ( a  -  ol  I 

U  Wo  'J  C  3t,  [XV[CUB  'J\ 


s  and 


[3D  3JJ  1  f-i  1  -  l Jjt  .,1  £ 
.Tc.  -iCvJ  [  eOlj  61[zd  J  3C,1 


(IV. 13) 


(IV. 14) 


(IV. 15) 


(IV. 16) 


(IV. 1?) 


(IV. 18) 


(IV. 19) 


(Equation  (IV. 19)  continued  on  next  page) 
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.  H  ^0  ♦  H  f  DD  *  1  "S* 

reo1  ^cv*  ro'L^CvJ  2t)c 


Everything  here  has  been  previously  obtained  except  for 


*CVX 


This  is  found  to  be 


VH  s  -JL  fo  f  <>£>  ^  v  J.  ^5 
'SC.,*  <<C^L  J  sTc  £CV  7>CV 


+  >/?  ole 

'SCv1- 


(IV. 20) 


In  this  Appendix,  we  have  obtained  the  extra  derivatives 

n 

.f.'IV'w  nrtl/Mil  n+i  ysrr  ^  ''  T  W  1V> 


necessary  for  calculating 


In  general  they  can  mostly 


be  expressed  in  terms  of  quantities  already  evaluated  in 


CONGRATS  for  ^  .  So  the  proper  combinations  of  the 

appropriate  quantities  were  programmed  into  CONGRATS,  along 

;\t© 

with  the  necessary  new  terms,  in  order  to  evaluate  x 

oCv 

This  was  done  by  Jean  Goertner  of  N'OL,  Once  this  was  accomplished, 


the  main  part  of  the  program  was  modified  to  calculate  the  pressure 
on  the  caustic  according  to  Equation  (4.1). 
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Appendix  V 


Input  Sound  Velocity  Data 


Profile  I  (Source  Depth  for  All  Calculations j  zQ 


Depth  Velocity 

m  m/s 


0 

1524.3 

457.2 

1521.57 

685.8 

1509.68 

853.44 

1495.35 

1036.3 

1490.48 

1188.7 

1490.46 

1478.28 

1491.7 

1981.2 

1499.01 

3352.8 

1516.99 

5250 

1548.78 

305  m) 


Profile  II  (  Source  Depth  for  All  Calculations!  zQ  =  305  m) 


Depth  Velocity 

m  m/s 


0 

1524.3 

457.2 

1521.57 

685.8 

1509.68 

853-44 

1495.35 

1036.3 

1490.48 

1188.7 

1490.46 

1478.28 

1491.7 

1981.2 

1499.01 

3352.8 

1516.99 

4572 

1537.42 

Profile  III  (Source  Depth  for  All  Calculations!  zQ  =  1000  m) 


Depth 

m 


Velocity 

m/s 


0  1570 
2250  1500 
5750  1570 


For  all  profiles,  p  =  p^  =  1.  Unit  consistency  is  not  necessary 
as  long  as  Equation  (5.12)  is  used  since  density  units  cancel. 
However,  for  evaluation  of  pressure  using  Equation  (5.9),  densities 

p 

should  be  in  MKS  units  in  order  that  pressure  be  in  nm  . 
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